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ABSTRACT 
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divergence  formulation,  and  a  semi-implicit  time  discretization.  In  the  first 
experiment  an  east-west  ridge  or  valley  is  placed  in  a  channel  with  east- 
west  periodic  conditions.  Linear  quasi-geostrophic  solutions  are  derived 
with  the  rigid  lid  assumption.  The  Rossby  waves  are  successfully  simulated 
in  the  model  with  linear  solutions  as  the  initial  conditions.  The  model  phase 
speeds  are  very  close  to  the  analytic  values  when  the  latter  are  properly 
corrected.  In  the  second  experiment  a  ridge  is  placed  across  the  channel  and 
the  Coriolis  parameter  is  set  to  zero.  The  initial  conditions  consist  of  a 
uniform  flow  through  the  channel  and  constant  free-surface  height.  The 
numerical  simulations  agree  with  hydraulic  jump  theory.  In  the  jump  cases 
the  model  predicts  increasing  wind  speeds  and  decreasing  free  surface 
heights.  Higher  spatial  resolution  would  be  required  to  properly  simulate 
the  details  of  the  hydraulic  jump  formation. 
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I.   INTRODUCTION 

A.     GENERAL 

Wilhelm  Bjerkens  (1904)  was  the  first  to  point  out  that  the  future 
meteorological  conditions  can  in  principle  be  obtained  by  an  integration  of 
differential  equations  which  govern  the  behavior  of  the  atmosphere.  Such 
an  integration  performed  using  numerical  methods  is  called  numerical 
^_S.3.Lhemr_s.rg_dmi£_Ho_n!_ 

Richardson  was  the  first  to  attempt  a  numerical  weather  prediction. 
After  very  long  and  time-consuming  computations,  he  obtained  a  totally 
unacceptable  result  (Richardson,  1922).  That  wrong  result,  and 
Richardson's  estimate  that  64,000  men  are  required  to  advance  the 
calculations  as  fast  as  the  weather  itself  is  advancing,  left  some  doubt 
about  the  practical  use  of  the  method.  However,  a  number  of  developments 
that  followed  improved  the  situation.  Mainly  due  to  the  work  of  Rossby  in 
the  late  1930' s,  it  became  clear  that  even  a  rather  simple  equation,  one  that 
describes  the  conservation  of  absolute  vorticity  following  the  motion  of  air 
particles,  suffices  for  an  approximate  description  of  large  scale  motions  of 
the  atmosphere. 

Finally,  in  1945,  the  first  electronic  computer,  ENIAC,  was 
constructed.  The  absolute  vorticity  conservation  equation,  and  ENIAC, 
were  used  by  Charney,  Fjortoft  and  Von  Neumann  in  the  late  1940' s  for  the 
first  successful   numerical    forecast    (Charney    et  al.    1950).    Much    faster 


computers,  and  improved  understanding  of  computational  problems,  now 
also  enable  long-term  integrations  of  the  basic  primitive  equations.  With 
the  introduction  of  each  new  generation  of  computers,  the  gap  between 
numerical  forecasts  and  atmospheric  observations  has  decreased,  but  the 
rate  at  which  this  gap  is  decreasing  appears  to  be  leveling  off.  This 
indicates  that  technological  improvements  in  computing  power  may  not  be 
the  primary  limitation  to  better  numerical  forecasts. 

During  the  past  15  years,  there  has  been  a  significant  effort  within  the 
numerical  weather  prediction  area  in  developing  limited-area,  fine-mesh 
primitive  equations  models  and  applying  them  to  operational,  short  range 
weather  forecasts.  An  important  practical  motivation  for  the  development  of 
regional  models  has  been  the  limited  success  of  operational  global  models 
in  the  prediction  of  precipitation  and  severe  weather.  A  parallel  motivation 
for  the  development  of  regional  models  is  their  potential  scientific  value  to 
researchers  studying  the  structure  and  dynamics  of  mesoscale  phenomena. 
(Keyser  and  Uccellini,  1987) 

B.     BACKGROUND 

There  are  two  fundamental  methods  of  simulating  the  atmospheric 
flow;  physical-models  and  mathematical-models  techniques.  With  the 
physical-models  technique,  we  construct  scale  model  replicas  of  observed 
ground  surface  characteristics  and  insert  them  into  a  chamber  such  as  a 
wind  tunnel.   The  flow   of  air  in   this   chamber   is   adjusted   so   as   to   best 


represent  the  large  scale,  observed  atmospheric  conditions.  Mathematical 
modeling,  on  the  other  hand,  makes  use  of  such  basic  analysis  techniques 
as  algebra  and  calculus  to  solve  directly  the  equations  governing 
atmospheric  flows. 

Numerical  solution  of  the  equations  of  motion  is  performed  using  the 
grid  point  method  for  most  applications.  Following  this  method  a  set  of 
points  is  introduced  in  the  region  of  interest  and  dependent  variables  are 
initially  defined  and  subsequently  computed  at  these  points.  This  set  of 
points  is  called  a  grid.  It  is  very  important  to  note  that  most  of  the  time  the 
grid  points  are  at  fixed  locations  in  the  horizontal.  This  means  that, 
according  to  the  Eulerian  system  of  equations,  space  and  time  coordinates 
are  chosen  as  independent  variables. 

With  the  grid  point  method,  the  most  common  way  of  solving  the 
governing  equations  is  to  find  approximate  expressions  for  derivatives 
appearing  in  the  equations.  The  required  approximate  expressions  are 
defined  using  values  of  the  dependent  variables  only  defined  at  the  grid 
points,  and  at  discrete  time  intervals.  Thus,  they  are  formed  using 
differences  of  dependent  variables  over  finite  space  and  time  intervals  ;  that 
is  the  reason  why  this  approach  is  called  the  finite  difference  method.  The 
approximations  for  derivatives  are  then  used  to  construct  a  system  of 
algebraic  equations  that  approximates  the  governing  partial  differential 
equations.  This  set  of  algebraic  equations  is  to  be  solved,  usually  using  an 
electronic  computer,  by  a  proper  step-wise  procedure  in  time. 


A  major  limiting  factor  of  finite  difference  approximations  is  the 
truncation  error.  That  is,  the  smaller  the  grid  interval,  the  smaller  the 
truncation  error.  For  a  finite  difference  model  to  improve  its  accuracy,  it 
would  require  increasing  the  grid  matrix  density.  This  would  require 
increased  computer  storage  and  computational  time.  The  problem  here  is 
not  a  simple  one.  It  goes  beyond  numerical  techniques  and  computer 
technology.  For  instance  if  we  further  reduce  the  grid  spacing  on  the  7LPE 
(National  Weather  Service  7  Layer  Primitive  Equation  Model)  we  do  not  get 
any  significant  improvement  in  the  accuracy  of  the  solution  (Woodward, 
1981).  The  required  additional  computer  capability  can  not  be  utilized  using 
finite  difference  methods.  Therefore,  new  numerical  integration  techniques 
must  be  investigated. 

Two  alternative  techniques,  the  spectral  method  and  the  finite  element 
method,  have  been  the  subject  of  intensive  research.  Both  the  spectral  and 
finite  element  methods  require  more  computational  time  per  forecast  time 
step  than  does  the  finite  difference  method.  That  is  why  both  the  spectral 
and  finite  element  methods  must  utilize  efficient  numerical  techniques  to  be 
considered  a  viable  option  for  numerical  weather  prediction.  For  long  range 
weather  predictions,  the  spectral  method  appears  to  be  a  natural  method, 
when  applied  over  the  globe  or  hemisphere.  This  is  due  to  the  existence  of 
efficient  transforms  for  the  nonlinear  terms  in  spherical  geometry. 
However,  the  spherical  harmonics  are  globally  rather  than  locally  defined. 


That  is   the  reason    why    the   finite   element   method   appears    to    be   more 
suitable  for  problems  of  more  detailed  limited  area  forecasting. 

The  Galerkin  Finite  Element  Method  (GFEM)  has  the  potential  to 
increase  efficiently  the  spatial  resolution  for  the  purpose  of  simulating 
accurately  the  small-scale  processes.  More  specifically,  the  GFEM  model 
used  by  Hinsman,  1983,  has  demonstrated  desirable  characteristics.  It 
propagates  atmospheric  waves  better  than  an  equivalent  finite  difference 
model.  It  also  allows  variable-resolution  grids  and  responds  better  than  an 
equivalent  finite  difference  model  near  the  smallest  gridlength.  Moving 
grids  can  be  achieved  with  no  apparent  noise  generation.  Finally,  it  can 
utilize  direct  solvers  and  is  a  natural  choice  for  vectorization  on  large 
computers.  Furthermore,  the  superior  small-scale  response  of  the  GFEM 
indicates  potential  increase  in  skill  for  regional  forecasting,  and  it  truly  is  a 
viable  option  for  simulation  of  atmospheric  flow. 

C.     OBJECTIVES 

The  main  objective  of  this  study  is  to  test  the  Galerkin  Finite  Element 
model  which  was  developed  by  Hinsman  (1983),  with  topography.  This 
shallow-water  model  uses  bilinear  basis  functions.  An  earlier  study  with 
the  same  model  was  carried  out  by  Neta  et  al.  (1986),  in  which  the  bottom 
topography  changes  linearly  to  the  north. 

In  this  thesis  two  types  of  bottom  topography  will  be  considered  in  a 
channel  domain   encompassed   by   solid   north-south   walls    with    east-west 


boundary  conditions.  In  the  first  case  the  bottom  is  composed  of  two 
regions.  These  regions  have  a  constant  but  opposite  northward  slope  so 
that  they  can  form  either  an  east-west  oriented  ridge  or  an  east-west 
oriented  valley  in  the  surface  topography.  Lines  of  constant  bottom  height 
run  parallel  to  the  x-axis,  and  pure  geostrophic  motion  is  possible  only  if 
the  v-component  of  the  velocity  is  identically  zero.  If  we  consider  no  mean 
flow  and  no  beta  effect,  topographic  Rossby  waves  can  exist  moving  either 
east  or  west  depending  on  the  slope  of  the  bottom. 

In  the  second  case,  a  mean  zonal  flow  passes  over  a  ridge  which 
extends  north-south  across  the  channel.  It  is  clear  now  that  the  lines  of 
constant  bottom  height  run  parallel  to  the  y-axis.  The  Coriolis  parameter 
for  this  case  is  set  to  zero,  and  the  formation  of  hydraulic  jumps  is  to  be 
investigated  for  different  values  of  the  mean  flow  and  peak  of  the 
topographic  ridge. 


IL.-THE  FINITE  ELEMENT  ME1HQD 

A.     HISTORICAL  INTRODUCTION 

The  finite  element  concept  is  to  a  large  extent  physical  rather  than 
abstract  in  nature,  and  it  has  been  used  in  a  variety  of  forms  for  centuries. 
The  basic  idea  has  always  been  to  replace  an  actual  problem  by  a  simpler 
one.  In  other  words  the  finite  element  method  is  the  replacement  of 
continuous  functions  by  piecewise  approximations,  usually  polynomials. 
Indeed  the  early  geometers  used  "finite  elements"  to  determine  an 
approximate  value  of  tc.  They  did  this  by  bounding  a  quadrant  of  a  circle 
with  inscribed  and  circumscribed  polygons,  the  straight-line  segments 
being  the  finite  element  approximations  to  an  arc  of  the  circle.  In  this  way 
they  were  able  to  obtain  extremely  accurate  estimates.  Upper  and  lower 
bounds  were  obtained,  and  by  taking  an  increasing  number  of  elements, 
monotonic  convergence  to  the  exact  solution  would  be  expected. 
Archimedes  used  these  ideas  to  determine  areas  of  plane  figures  and 
volumes  of  solids,  although  of  course  he  did  not  have  a  precise  concept  of 
a  limiting  procedure.  The  interesting  point  here  is  that  while  many 
problems  of  applied  mathematics  are  posed  in  terms  of  differential 
equations,  the  finite  element  solution  of  such  equations  utilizes  ideas  which 
are  in  fact  much  older  than  those  used  to  set  up  the  equations  initially. 

The    modern    use   of    finite    elements    really    started    in    the    field    of 
structural    engineering.    Probably    the    first    attempts    were    by    Hrennikoff 


(1941)  and  McHenry  (1943)  who  developed  analogies  between  actual 
discrete  elements  and  the  corresponding  portions  of  a  continuous  solid.  The 
term  "finite  element"  was  introduced  later  by  Clough  (1960)  in  a  paper 
describing  an  application  in  plane  elasticity. 

The  engineers  had  put  the  finite  element  method  on  the  map  as  a 
practical  technique  for  solving  their  elasticity  problems,  and  although  a 
rigorous  mathematical  basis  had  not  been  developed,  the  next  few  years 
saw  an  expansion  of  the  method  to  solve  a  large  variety  of  structural 
problems.  The  workers  in  the  early-1960s  soon  turned  their  attention 
towards  the  solution  of  non-linear  problems.  Turner  et  al.  (1960)  showed 
how  to  use  an  incremental  technique  to  solve  geometrically  non-linear 
problems,  i.  e.,  problems  in  which  the  strains  remains  small  but  the 
displacements  are  large.  Stability  analysis  also  comes  into  this  category  and 
was  discussed  by  Martin  (1965).  Plasticity  problems,  involving  non-linear 
material  behavior  were  modelled  at  this  time  (Gallagher  et  al.  1962),  and 
the  method  was  also  applied  to  the  solution  of  problems  in  visco-elasticity 
(Zienkiewicz  et  al.  1968). 

Finally,  besides  the  static  analysis,  dynamic  problems  were  also  being 
tackled,  and  Archer  (1963)  introduced  the  concept  of  the  consistent  mass 
matrix.  Both  vibration  problems  (Zienkiewicz  et  al.  1966)  and  transient 
problems  (Koenig  and  Davids,  1969)  were  considered.  Thus  the  period 
from  its  conception  in  early-1950s  to  the  mid-1960s  saw  the  method  being 
applied  extensively  by  the  engineering   community.    Once   it   was    realized 


that  the  method  could  be  interpreted  in  terms  of  variational  methods,  the 
mathematicians  and  engineers  were  brought  together,  and  many  extensions 
of  the  method  to  new  areas  soon  followed.  In  particular  it  was  realized  that 
the  concept  of  piecewise  polynomial  approximation  offered  a  simple  and 
efficient  procedure  for  the  application  of  the  classical  Rayleigh-Ritz 
method.  The  method  had  now  become  an  important  technique  from  both  a 
practical  and  theoretical  point  of  view,  and  the  number  of  published  paper 
using  this  method  began  to  increase  at  a  tremendous  rate. 

In  the  area  of  meteorology  also,  the  finite-element  method  has  been 
successfully  employed  in  the  horizontal  representation  of  atmospheric 
variables  in  numerical  weather  prediction  and  atmospheric  modeling 
(Cullen,  1974a,  1974b,  1979;  Hinsman,  1975,  1983;  Staniforth  and 
Mitchell,  1977,  1978).  The  finite  element  method  when  applied  to 
meteorological  equations  gives  very  accurate  phase  propagation  and  also 
handles  nonlinearities  very  well. 

B.     FINITE  ELEMENT  APPROXIMATIONS 

In  contrast  to  the  finite  difference  schemes  wherein  the  domain  of 
interest  is  replaced  by  a  set  of  discrete  points,  in  the  finite  element  method 
the  domain  is  divided  into  subdomains  called  finite  elements. The  unknown 
function,  let  us  name  it  u,  is  represented  within  each  element  by  an 
interpolating  polynomial  which  is  continuous  along  with  its  derivatives  to  a 
specified  order  within  the  element.  Generally,  the  interpolating  function  is 


of  lower-order  continuity  between  elements  than  within  an  element.  Thus 
the  fundamental  building  block  in  the  finite  element  method  is  the 
subdomain  or  element. 

C.     THE  METHOD  OF  WEIGHTED  RESIDUALS 

There  are  several  ways  that  can  lead  us  to  the  same  finite  element 
formulation.  A  conceptually  simple  approach  can  be  formulated  using  the 
method  of  weighted  residuals.  The  two  primarily  special  cases  of  the 
method  of  weighted  residuals  (MWR)  are  the  Galerkin  and  the  collocation 
methods. 

In  the  method  of  weighted  residuals,  the  desired  function  u  (•)  is 
replaced  by  a  finite  series  approximation  as 

N 

u  (•)  =  u  (•)  =  ]£  U.  <J>.  (•).  (  2  -  1 ) 

j  =  i 

In  general,  the  set  of  functions  <J>j  (•),  j  =  1,2, ...,N,  can  be  defined  over 
both  the  time  and  space  domain  and  Uj,  j  =  1,2. ..,N,  are  undetermined 
coefficients.  The  equation  (2-1)  can  be  written 

N  -  1 

u(.)  =  U0(|)0(-)  +  XUJ(I)J(,)'  (2"2) 

i-i 
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where  <|>j  (•)  satisfy  the  homogeneous  boundary  conditions.  The  functions 
<j)j  (•)  are  chosen  to  be  polynomials  that  satisfy  certain  of  the  boundary 
conditions  imposed  on  the  problem.  These  functions  are  variously  denoted 
shape  functions,  basis  functions,  and  interpolation  functions,  depending 
upon  the  discipline  in  which  the  method  is  being  applied.  Even  if  we 
choose  the  basis  functions  to  satisfy  all  boundary  conditions,  they  will  not 
normally  satisfy  the  PDE  as  well.  If  we  now  substitute  the  u  (•)  into  the 
PDE,  say  Lu  -  f  =  0,  we  can  easily  get 

Lu(-)-f  =  R(.),  (2-3) 

where  R  (•)  is  a  residual. 

Our  objective  at  this  point  should  be  to  select  the  undetermined 
coefficients  Uj  such  that  this  residual  is  minimized  in  some  sense.  A 
straightforward  scheme  would  be  to  set  the  integral  of  R  (•)    to  zero  as 

J  j  R  (•)  dv  dt  =  0.  (  2  -  4  ) 

t      V 

This  scheme,  however,  generates  only  one  equation  for  the  N  unknown 
coefficients  Uj.  It  can  be  suitably  modified  by  introducing  weighting 
functions  wj  (•),  i  =  1,2,. ..,N. 

Setting  the  integral  of  each  weighted  residual  to  zero  yields  n 
independent  equations 
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J  J  R  (•)  w.  (•)  dv  dt  =  0,  i  =  1,2,...,N.  (  2  -  5  ) 

t.       V 


At  this  point  we  can  solve  equation  (2-5),  in  theory  at  least,  for  N 
unknown  coefficients.  Equation  (2-5)  represents  the  general  equation 
describing  the  MWR,  and  a  multiplicity  of  schemes  arise  out  of  this  one 
expression  through  the  definition  of  the  weighting  functions  wj. 

Among  the  MWR  family  of  methods,  the  Galerkin,  subdomain,  and 
collocation  schemes  are  most  commonly  encountered  in  practice.  The  one- 
dimensional  weighting  function  for  the  Galerkin  scheme  is  illustrated  in 
Fig.  2.1,  for  the  subdomain  scheme  in  Fig.  2.2,  and  for  the  collocation 
scheme  in  Fig.  2.  3. 

D.     GALERKIN   METHOD 

The  Galerkin  method  results  when  the  weighting  function  is  chosen  to 
be  the  basis  function,  as  defined  in  (2-1).  Thus  we  have 

J  J  R  (•)  <».  (•)  dv  dt  =  0,  i  =  1,2,...,N.  (  2  -  6  ) 

t       V 

The  basis  functions  are  formally  required  to  be  members  of  a  complete  set 
of  functions.    Because   a  complete   set  of  functions   can   exactly   represent 
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Fig.  2.1  Schematic  representation  of  the  one-dimensional 
weighting  functions  for  the  Galerkin  method.  (It  is 
assumed  that  the  chapeau  function  is  used  as  a  basis 
function). 
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Fig.   2.  2    As  in  Fig.  2. 1,  except  for  the  subdomain  method. 
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Fig.   2.  3    As  in  Fig.  2. 1,  except  for  the  collocation  method. 
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any  function  of  a  given  class,  the  series  of  (2-1)  is  inherently  capable  of 
representing  the  exact  solution  as  the  number  of  terms  in  the  series  is 
increased. 

The  requirement  of  completeness  allows  an  alternative  interpretation  of 
the  Galerkin  formulation.  A  continuous  function  must  be  zero  if  it  is 
orthogonal  to  every  member  of  a  complete  set.  The  Galerkin  method  can  be 
viewed  as  a  scheme  in  which  the  residual  is  forced  to  zero  in  the  sense  that 
it  is  made  orthogonal  to  the  complete  set  of  functions  (^  (•). 

In  Fig.  2. 1  the  weighting  function  (and  therefore  the  basis  function)  is 
a  hatshaped,  piecewise  linear  function,  and  because  of  its  hatlike 
appearance,  it  is  sometimes  called  a  "chapeau"  function.  It  is  often 
encountered  in  the  formulation  of  the  finite  element  method.  The  chapeau 
function  is  the  simplest  of  the  basis  functions  in  common  use. 

E.     TWO-DIMENSIONAL  BASIS   FUNCTIONS 

The  extension  of  the  weighted  residual  method  to  higher  dimensions  is 
relatively  straightforward,  provided  that  regular  rectangular  subspaces  are 
employed.  We  can  easily  visualize  the  two-dimensional  case  using  an 
elementary  example,  (Lapidus  and  Pinder,  1982).  The  discretized  domain  is 
simply  illustrated  in  Fig.  2.4,  and  the  two-dimensional  basis  function  that 
is  linear  along  each  side  in  Fig.  2.5. 
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Fig.  2.4  Discretized  domain  for  two-dimensional  heat  flow. 
Finite  element  nodes  are  indicated  by  small  circles,  and 
the  elements  themselves  by  Roman  numerals. 
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Fig.    2.5    Two-dimensional    basis    function    that    is    linear    along 
each  side. 
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Consider  now  the  problem  of  two-dimensional,  time-independent  heat 
flow  in  a  rectangular  plate  with  a  heat  source  located  at  the  center  of  the 
plate.  The  governing  equations  for  this  problem  will  be  given  by 


£(T)  =  T      +  T     =Q  (2-7) 

v       '  xx  yy       ^  v  ' 


T  (x,2)  =  1  (  2  -  8  ) 


T  (0,y)  =  1  (  2  -  9  ) 


Ty  (x,0)  =  0  (  2  -  10  ) 


Tx(2,y)  =  0  (2-11) 

Q(x,y)  =  Qw8(x-l)5(y-l),  (2-12) 

where  x  and  y  are  Cartesian  coordinates,  Qw  is  the  heat  source,  and  5  is  the 
Dirac  delta  function. 

Let  us  first  define  the  trial  functions: 


T  (x,y)  -  f  (x,y)  =  ]£  T.  <)>.  (x,y).  (  2  -  13  ) 

j=  i 
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For  convenience,  it  is  advantageous  to  define  basis  functions  in 
local  (^,tj)  coordinates  where  upon  (2-13)  becomes,  for  £QQ.h_£lement± 

4 

T  (x,y)  -  f  (x,y)  =  £  T.  <j>.  <&J\)  (  2  -  14  ) 

j-i 

and  <J»j  (^,TJ)  are  bilinear  chapeau  functions  such  as  those  illustrated  in  Fig. 

2.5. 

2*_JLl£iLJI 

At  this  point  we  can  easily  formulate  the  integral  equations  using 
Galerkin's  procedure.  We  can  visualize  the  whole  procedure  as  the 
requirement  of  orthogonality  between  the  residual  R  and  each  basis 
function,  as 

J  R  (x,y)  4>.  (x,y)  dx  dy  =  0,  i  =  1,...,9,  (  2  -  15  ) 

A 


where 


R(x,y)s£(T)-Q.  (2-16) 

A  is  the  domain  of  integration  and  £  is  defined  by  (2-7).  Substituting  (2-7), 
and  (2-16)  into  (2-15)  yields 
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J  <  Txx  +  Tyy  "  Q  )  *i  (x'y}  &  ^  =  °'  *  =  1""'9-  (  2  -  17  ) 


Applying  Green's  theorem  to  (2-17)  to  incorporate  second-  and  third-type 
boundary  conditions  directly  into  the  set  of  integral  equations,  equation  (2- 

17)  becomes 

f(T6.  +  T(b.  +  Q(|).  )  dx  dy  -  f  ( f    1  +  f    1   )  <f>.  ds  =  0, 

j  v    x  Txi       y  Yyi     x  Yi '  J       xx        yy1 

A  S 

i=l,...,9,  (2-18) 

where  lx  and  ly  are  direction  cosines  with  respect  to  the  normal  to  the  curve 
S,  the  boundary  of  the  domain  A.  In  this  problem,  the  second  term  of  (2- 

18)  will    be    used    to    conveniently    define    the    zero-gradient    Neumann 
boundary  conditions  of  (2-10)  and  (2-11). 

If  we  now  substitute  the  trial  functions,  as  defined  by  (2-13)  into 
(2-18)  we  can  obtain  the  following  set  of  algebraic  equations 


J(  JvMxi  +  ^V  +  Q(,,i)dxdy 


A      J-l 


-J  (  f  x  lx  +  f  y  ly  )  <>.  ds  =  0,  i  =  1,...,9.  (  2  -  19  ) 
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Because  <j>j  is  defined  such  that  it  is  nonzero  only  over  elements  adjacent  to 
mode  i  (see  Fig.  2.4  and  2.5)  the  integrations  of  (2-19)  may  be  performed 
piecewise  over  each  element  and  subsequently  summed.  Thus  we  can  write 


4  ,    9 


X   J(I,V^  +  ^)  +  QOdxdy 


e-l       a         J=l 


-  J  (f  x  lx  +  f  y  ly  )  ^  ds  =  0,  i  =  1.....9,  (  2  -  20 ) 


se 


where  Ae  is  the  area  of  element  e,  and  Se  is  the  curve  bounding  Ae. 

It  is  now  time  to  formulate  the  matrix  equation.  In  general,  the 
best  way  is  to  express  (2-20)  in  terms  of  the  local  (£,,T[)  coordinate  system 
to  facilitate  integration  (see  Fig.  2.6a,  Fig.  2.6b).  This  is  easily 
accomplished  provided  that  the  relationship  between  derivatives  of  <j>j  in 
each  coordinate  system  is  readily  available.  To  find  this  relationship,  we 
have  to  employ  the  chain  rule  to  obtain 
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Fig.    2.6a  Rectangular  element  in  (x,y)  coordinates. 
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Fig.    2.6b  Same  element  as    Fig.    2.6a  transformed   into   the 

(£,,T\)  local  coordinate  system. 
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^j  _<tyj_3x      jfoj  dy 
"^""9x  a^      dy  a^  ' 


This  set  of  equations  can  be  written  as 


L(Vn. 


(x)        (y) 

"H  11. 


(0-) 


VYj/y 


or 


J    X 


(40 

.    J  yJ 


=[J] 


-l 


(2-21) 


where  [J]  is  the  Jacobian  matrix.  For  our  problem  we  can  easily  evaluate 
the  [J]  as 


M- 


0.5      0.0 
0.0      0.5 


and  the  [J]-1  as 
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pr- 


2.0      0.0 
0.0      2.0 


(  2  -  22  ) 


Equations  (2-20)  and  (2-22)  may  be  combined,  which  can  give  us 


4     +}  +}  9 


£J  j(£v4v*+4w+Q*iKd^dTi 

e=l.i    -1        j=  1 

-  I  (f    1  +f    1  )<J».ds  =  0,  i=  1,...,9. 


(  2  -  23  ) 


In     changing     the     limits     of     integration     we     introduce     the     following 
relationship 


dx  dy  =  det  [J]  d^  drj. 


(2-24) 


Because  in  our  problem  the  elements  are  all  of  the  same  size,  the 
four  integrals  appearing  in  each  element  in  (2-23)  are  identical  for  each 
element.  The  assembly  process,  that  is,  the  transformation  from  element  to 
global  integrations,  can  be  easily  visualized  using  Table  I  which  provides 
us  with  the  relationship  between  Local  Node  Numbers  and  Global  Node 
Numbers. 

The  assembly  procedure  now  calls  for  extracting  information 
concerning    the   integration    at   the   element   level   from   the    Table    I.    From 
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Element  Nodes       ""*■"" — -«^iiii 

Global  Nodes 
(Global  Node  Numbers) 

(Local  Node  Numbers) 

I 

II 

III 

IV 

1 
2 

3 
4 

1 
4 
5 

2 

2 
5 
6 
3 

4 
7 
8 
5 

4 
8 
9 
6 

Table  I       Relationship  between  Local  Node  Numbers  and  Global 
Node  Numbers. 
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column  I  we  can  see  that  global  nodes  4  and  5  correspond  to  element  nodes 
2  and  3  in  element  I.  In  order  to  find  the  contribution  to  the  global  integral 
concerning  these  nodes,  we  have  to  create  Table  II  which  evaluates  the 
integrals  for  one  element  in  local  coordinates  using  equation  (2-23).  From 
Table  II  the  seeking  contribution  is  found  in  row  2,  column  3,  and  it  is 
equal  to  -1/6.  This  value  should  be  placed  in  matrix  location  (4,5)  of  the 
global  matrix.  However,  this  value  is  not  final,  because  there  is  also 
information  to  be  retrieved  from  element  III.  This  information  is  located  in 
the  position  (1,4)  of  the  element  matrix  of  Table  II,  and  it  is  equals  to  -1/6. 
This  value  has  to  be  summed  with  the  previous  integral  value  and  the  new 
value  should  to  be  placed  in  the  same  global  matrix  position  (4,5). 
Combination  of  the  two  integral  values  yields  the  final  value  of  -1/3,  which 
can  be  seen  in  Table  III  in  row  4,  column  5. 

In  general,  the  element  coefficient  matrix  (Table  III)  is  different 
for  each  element  because  of  changes  in  either  element  geometry  or 
parameter  values.  Moreover,  it  is  often  necessary  to  perform  the 
integrations  of  (2-23)  numerically. 
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X 

1 

2 

3 

4 

1 

2/3 

-1/6 

-1/3 

-1/6 

2 

-1/6 

2/3 

-1/6 

-1/3 

3 

-1/3 

-1/6 

2/3 

-1/6 

4 

-1/6 

-1/3 

-1/6 

2/3 

Table  II     Evaluated      Integrals      for      One      Element      in      Local 
Coordinates  for  equation  (2-20). 
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\ 

1 

2 

3 

4 

5 

6 

7 

8 

9 

1 

2/3 

-1/6 

0.0 

-1/6 

-1/3 

0.0 

0.0 

0.0 

0.0 

2 

-1/6 

4/3 

-1/6 

-1/3 

-1/3 

-1/3 

0.0 

0.0 

0.0 

3 

0.0 

-1/6 

2/3 

0.0 

-1/3 

-1/6 

0.0 

0.0 

0.0 

4 

-1/6 

-1/3 

0.0 

4/3 

-1/3 

0.0 

-1/6 

-1/3 

0.0 

5 

-1/3 

-1/3 

-1/3 

-1/3 

8/3 

-1/3 

-1/3 

-1/3 

-1/3 

6 

0.0 

-1/3 

-1/6 

0.0 

-1/3 

4/3 

0.0 

-1/3 

-1/6 

7 

0.0 

0.0 

0.0 

-1/6 

-1/3 

0.0 

2/3 

-1/6 

0.0 

8 

0.0 

0.0 

0.0 

-1/3 

-1/3 

-1/3 

-1/6 

4/3 

-1/6 

9 

0.0 

0.0 

0.0 

0.0 

-1/3 

-1/6 

0.0 

-1/6 

2/3 

Table   III  Evaluated  Integrals  for  the  Four  Elements  of  Fig.    2.4 
in  Global  Coordinates  for  Equation  (2-20). 
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We  now  solve  the  matrix  equation  (2-23)  using  Table  III  as 


2/3     -1/6     0.0     -1/6     -1/3     0.0     0.0     0.0     0.0 


1/6      4/3    -1/6     -1/3     -1/3    -1/3     0.0     0.0     0.0 


0.0     -1/6     2/3      0.0     -1/3    -1/6     0.0     0.0     0.0 


1/6     -1/3     0.0      4/3     -1/3     0.0    -1/6   -1/3     0.0 


-1/3     -1/3    -1/3     -1/3      8/3    -1/3    -1/3    -1/3    -1/3 


0.0     -1/3    -1/6      0.0     -1/3     4/3     0.0    -1/3    -1/6 


0.0      0.0     0.0     -1/6     -1/3     0.0     2/3    -1/6     0.0 


0.0      0.0     0.0     -1/3     -1/3    -1/3    -1/6     4/3    -1/6 


0.0      0.0     0.0      0.0     -1/3    -1/6     0.0   -1/6     2/3 


T2 


T3 


T4 


T5 


T, 


T7 


T9 
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I  ( f    1  +  f    1   )  0.  ds 

s 

I  ( f    1  +  f    1  )  A,  ds 

J    v       xx  yy/T2 

S 

f(f     1    +f     1    )(j),  ds 

jv     xx         y  y      3 


0.0 


-Qw 

f(f    1  +f    1  )<b.  ds 

jv     xx         yyo 


0.0 


0.0 


f(t    1  +f    1  )(|>Qds 

J    v       xx  yy/T9 


(  2  -  25  ) 
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The  line  integrals  appearing  on  the  right-hand  side  of  (2-25) 
represent  flux  boundary  conditions.  In  general,  when  Dirichlet  boundaries 
are  specified,  it  is  necessary  to  expand  and  evaluate  the  line  integrals. 
However,  when  we  use  functions,  that  have  continuous  derivatives  up  to 
the  Oth  order  as  bases,  one  can  condense  rows  containing  the  known 
temperature  values  out  of  the  matrix  equation.  The  reduced  matrix  equation 
is 


where 


[A]  *  {b}  =  {f}. 


(  6  -  26  ) 


[A]. 


4/3 

-1/3 

-1/6 

-1/3 

-1/3 

8/3 

-1/3 

-1/3 

-1/6 

-1/3 

2/3 

-1/6 

-1/3 

-1/3 

-1/6 

4/3 

(  6  -  27  ) 


T. 


Tc 


{b}  = 


(  6  -  28  ) 
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{f}  = 


1/2 

0.0 
1/2 


(  6  -  29  ) 


Solving  (2-26)  directly  we  can  easily  get 


T4 

T5 

T7 

T8 


0.783 


0.525 


0.655 


0.783 


(  2  -  30  ) 


The  coefficients  T4,  T5,  T7,  and  Ts  represent  the  values  of  the 
temperature  at  nodes  4,  5,  7,  and  8  because  the  basis  functions  are  defined 
such  that  <J>i  is  unity  at  node  i  and  zero  elsewhere. 

At  the  same  time  the  analytical  solution  for  our  problem  is 


T  =  U  +  V  +  W, 


(2-31) 
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where 


K   n  =  0 


sin[(n  +  l/2)7c(a  -  y)/a]  cosh[(n  +  l/2)jc(a  -  x)/a] 
(n  +  1/2)  cosh(n  +  1/2)71 


(  2  -  32  ) 


«   n  =  0 


cosh[(n  +  l/2)juy  /  a]  sin[(n  +1/2)jcx  /  a] 
(n  +  1/2)  cosh(n  +  1/2)7C 


(  2  -  33  ) 


2  ^  sin[(n+l/2)7Cx/a]sin[(n+l/2)7C^/a]sinh[(n+l/2)7c(a-y)]cosh[(n+l/2)7cri/a] 


^n  =  0 


(n  +  1/2)  cosh[(n  +  l/2)7c] 


(2-34) 


where  a  is  the  length  of  the  side  of  the  square  (in  our  case  equals  to  2), 
and  £,  tj  are  the  x  and  y  locations  of  the  heat  source,  respectively.  The 
analytical  solution  yields 


"V 

"0.756" 

T2 

0.127 

T3 

0.719 

A 

0.756 

(  2  -  35  ) 
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It  is  obvious  that  the  numerical  solution  at  the  singular  point  (1,1),  which 
corresponds  to  the  nodal  location  5,  is  not  very  accurate.  This  is  not  odd 
since  we  are  attempting  to  represent  a  rapidly  varying  function  with  only 
four  bilinear  elements.  As  we  mentioned  before  the  series  of  (2-1)  is 
inherently  capable  of  representing  the  exact  solution  as  the  number  of  terms 
in  the  series  is  increased. 
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hil_the;  g  h allow- water  mqdel  and  the  topographic 

EHaS.BY_SAXE 

A.  GENERAL 

In  order  to  study  motions  of  atmospheric  and  oceanic  relevance,  we 
can  use  a  shallow,  rotating  layer  of  homogeneous  fluid  which  is 
incompressible,  and  inviscid.  This  model  (Shallow-Water  Model)  ignores 
completely  the  presence  of  stratification,  but  experience  has  shown  that  it 
is  capable  of  describing  important  aspects  of  atmospheric  and  oceanic 
motions.  Because  of  this,  it  is  useful  to  deal  with  shallow  fluid  systems  in 
trying  to  determine  general  principles  of  hydrodynamical  behavior  of  the 
atmosphere. 

The  major  physical  characteristics  concerning  the  shallow-water 
model,  are  found  to  apply  well  for  more  complex  systems.  If  we  were  to 
analyze  the  types  of  wave  motions  associated  with  more  complex  forms  of 
the  primitive  atmospheric  equations,  (e.g.,  with  vertical  stratification, 
compressibility,  etc.),  we  would  find  essentially  the  same  types,  namely 
gravity-inertia  waves  and  Rossby  waves,  and  deep  and  shallow  motions 
would  exist  simultaneously. 

B.  THE  SHALLOW-WATER  MODEL 

Let  us  consider  a  sheet  of  fluid  with  constant  and  uniform  density  as 
illustrated  in  Fig.   3.1.  The  height  of  the  surface  of  the  fluid  is   given  by 
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ft  f 

P--  CONSTANT 

1  _ 

g 

L                           -i 

1  -* 

/  y,v               h(x,y,t) 

/ 
/ 

1  hB  (x,y) 

/                                                                                                                                                 

_- — ■ 

fc- 

1 

K,U 

Fig.   3.1    The  Shallow-Water  Model. 
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h(x,y,t).  We  also  model  the  body  force  arising  from  the  potential  (J)  as  a 
vector,  g,  normal  to  the  z  =  0  surface.  The  z-axis  coincides  with  the 
rotation  axis  of  the  fluid,  so  that  in  this  particular  case  the  Coriolis 
parameter  f  is  simply  given  by  2Q..  The  rigid  bottom  topography  is  defined 
by  the  surface  z  =  hB(x,y),  which  is  not  a  function  of  time  (t).  Finally,  we 
assume  that  the  fluid  is  inviscid  (|i=0),  that  is,  only  motions  for  which 
viscosity  is  unimportant  are  considered. 

Also  we  suppose  that  a  characteristic  value  for  the  depth  can  be 
sensibly  chosen,  say  D,  and  D  also  characterizes  the  vertical  scale  of  the 
motion  as  well.  In  the  same  way  we  consider  that  a  characteristic  horizontal 
length  scale  for  the  motion  exists,  which  we  call  L.  The  fundamental 
relationship  which  characterizes  the  shallow-water  theory  is  given  by 


5  =  £  «  1.  (3-1) 


Here,  it  is  important  to  note  that  the  fluid  is  rotating,  so  that  Coriolis 
accelerations  can  be  important.  The  fluid  layer  is  flat  rather  than  forming  a 
spherical  shell,  and  its  major  physical  deficiency  as  mentioned  above  is  the 
absence  of  the  density  stratification  which  is  present  in  the  real 
atmosphere. 
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C.     SMALL-AMPLITUDE  MOTIONS 

Dealing  with  the  Shallow-Water  Model,  and  since  by  hypothesis  5<<1, 
we  are  able  to  express  the  hydrostatic  approximation  in  the  following  form: 

dp 

£-pg.  (3-2) 


We  can  integrate  equation  (3.2)  resulting  in 

p  =  -  p  g  z  +  A(x,  y,  t).  (  3  -  3  ) 

Applying  now  the  obvious  boundary  condition 

p(x,y,h)  =  p0,  (  3  -  4  ) 

where  po  is  a  constant,  we  can  easily  get 

p  =  pg(h-z)  +  p0.  (  3  -  5  ) 

At  this  stage,  we  can  clearly  see  that  the  horizontal  pressure  gradient 
is  independent  of  z,  i.e. 

/=pgy-,  (3-6) 

ox  ox 
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aTpg97-  (3"7) 


in  such  way  that  the  horizontal  accelerations  must  be  independent  of  z.  That 
is  why  we  can  also  assume  that  the  horizontal  velocities  themselves  remain 
independent  of  z,  if  they  are  so  initially.  Applying  now  the  Taylor- 
Proudman  theorem  we  are  able  to  write  the  horizontal  momentum  equation 
as 


3u        3u       3u     „  3h  ,  „    „ 

3T+u^r+varfv  =  -ga^  (3"8) 


9v        dv       dv      „  3h  ,  „    « , 

aT+uar+v87+fu  =  -g3F-  (3"9) 


The  specification  of  incompressibility  for  the  Shallow-Water  Model, 
decouples  the  dynamics  from  the  thermodynamics  and  reduces  the  equation 
of  mass  conservation  to  the  condition  of  incompressibility  given  in  the 
following  form: 

du      dv     9w 

a7+ar<fe=a  (3-10) 

Having  in  mind  that  u  and  v  are  independent  of  z,  equation   (3-10)  can  be 
integrated  in  z  as 
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,du      9vN 
w(x,y,t)  =  -  z(-^-  +  j-)  +  <5(x,y,t).  (3-11) 


Equation  (3-11)  can  be  further  manipulated  using  the  condition  of  no 
normal  flow  at  the  rigid  surface  z  =  hs  resulting  in 

g.  +  A((h-hB)u)+|r{(h-hB)v}=0.  (3-12) 


If  we  now  consider  the  total  depth  (H)  is  given  by 


H  =  h-hB,  (3-13) 


then  the  equation  of  mass  conservation  (3-12)  becomes 


— +  _(uH)  +  —  (vH)  =  0.  (3-15) 


Our  next  step  now  should  be  to  linearize  the  set  of  equations  (3-8),  (3- 
9),  and  (3-15),  by  studying  small-amplitude  motions.  This  is  very 
important  because  the  presence  of  solutions  representing  free  oscillations 
often  demonstrates  fundamental  mechanisms  which  occur  in  more 
complicated  situations  as  we  mentioned  before. 

Let  us  now  consider  the  thickness  of  the  fluid  layer  in  absence  of 
motion  to  be  given  as 
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H(x,y,t)  =  HQ(x,y)  +  Ti(x,y,t),  (  3  -  16  ) 

where  r\  «  Ho.  Furthermore 

»  u„  .  Vuu,  (  3  -  17  ) 


at        H  *    H' 

or  in  other  words  u  and  v  are  considered  small  enough.  Linearizing  now  the 
equations  (3-8),  (3-9),  and  (3-15),  we  can  obtain 


£-*-«£■ 


§>=-#■ 


^+a7(uHo)  +  a?(vHo)  =  0-  (3"20) 


where  all  the  quadratic  terms  in  the  dynamical  variables  u,  v,  rj  with 
respect  to  the  linear  terms  are  ignored.  If  we  now  define  the  linearized 
mass  flux  vector  given  by 

U  =  iU+jV,  (3-21) 


where 
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U  =  uHQ,  (3-22) 


V  =  v  HQ,  (  3  -  23  ) 


equations  (3-18),  (3-19),  and  (3-20)  become 


£-"-*£.  <>■*> 


IJUfU-gH^.  (3-25) 


^+87+a7=a  (3"26) 


At  this  point  it  is  possible  for  us  to  obtain  an  equation  in  the  single 
variable  r\  as 


^  [(\+  f2)  r\  -  V  .  (CqVti)]  -  g  f  J  (H0,ti)  =0,  (  3  -  27  ) 

at 


where 


Co=gH0,  (3-28) 
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JWnrw"sr£-  (3"29) 


It  can  be  seen  that  the  Jacobian  term  is  the  effect  of  the  geostrophic 
wind  blowing  across  isobaths.  For  low  frequency  motion  and  small 
variation  in  Ho,  the  first  term  in  (3-27)  represents  the  time  rate  of  change 
of  the  quasigeostrophic  vorticity.  The  motion  that  results  from  this  kind  of 
balance  is  a  topographic  Rossby  wave  which  we  examine  in  more  detail  in 
the  next  sections.  The  velocities  components  u  and  v  can  be  also  found  in 
terms  of  T|,  given  as 

D.     THE  TOPOGRAPHIC  ROSSBY  WAVE 

Following  Pedlosky  (1987),  let  us  consider  that  Ho  in  equation  (3-27) 
varies  slightly  in  the  y-direction    given  by 


HQ  =  D-Ds|-,  (3-32) 


where  D  is  a  constant, 
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s  is  the  slope  of  Ho  in  the  y-direction,  and 
L  is  the  width  of  the  channel. 
In  this  particular  case  (Fig.  3.2),  pure  geostrophic  motion  is  possible 
only  if  v  is  zero.   Also,   lines   of  constant  Ho   are  parallel   to   the  x-axis. 
Motion  across  the  isobaths  of  fluid  columns  will  cause  them  to  stretch  or 
contract.  Therefore,  there  is  a  possibility  of  a  different  mode  of  motion  to 
exist  depending  on  the  combined  effect  of  rotation  and  bottom  slope. 
Assuming 

s  «   1,  (3-33) 

and  r\  to  be  of  the  following  form 

ti  =  Re  {  lf(y)  exp[i(kx  -  at)]  }  ,  (  3  -  34 ) 

substitution  in  (3-27)  yields 


2  *  2-2 

y   d  ti      c  dTi         *  a  -  F      2  y      fs 

(1-sf)— L.i     '    +T,   [_^-k2(l-sf) k]=0,  (3-35) 

L    dy2      L  dy  gD  L      La 


with  the  following  boundary  conditions: 
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Fig.   3.2     The  infinite  channel  of  width  L,  rotating  with  angular 
velocity  f/2. 
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dn       f  k    * 

-T-  +  —  Tl   =0,  (3-36) 

dy     a 


on  y  =  0,  L. 

If  we  now  assume  that 


(l-s£)=l,  (3-37) 


which  appears  to  be  a  very  realistic  approximation,  then  (3-35)  becomes 

2    *  *  2    2 

d-n      sdr|       ra-i,2fs,.*-  ,  ~    -« * 

— L. -__J_+[_ — k  -- — km   =0.  (3-38) 

dy2     L  **         C20  L  a 

Here  only  in  terms  where  s  can  be  compared  with  quantities  of  order 
unity  has  it  been  neglected.  Solving  (3-38)  we  can  get 


Tl*  =  exp(-^-)  [A  sin(ay)  +  B  cos(ay)],  (  3  -  39  ) 


where  a  is  given  by 


/  &  - 1      „  2       s    „     fks  ,  „     ._ . 

a=      I— (k  +— -) .  (3-40) 

S  4L       cjL 
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Applying   now   the  boundary   conditions   given   by   (3-36)   we  can   get   the 
following  eigenvalue  problem 

(o2  -  f2)  (ct2  -  k2C^)  sin  (aL)  =  0.  (  3  -  41  ) 

It  is  obvious,  that  to  the  lowest  order  in  s,  the  slope  does  not  alter  the 
Kelvin  mode  (that  is,  because  the  factors  multiplying  sin  (aL)  are  the  same 
as  for  the  case  of  a  flat  bottom).  The  roots  corresponding  to  the  zeros  of 
sin  (aL)  are  given  by 

a2 »-C20(k2  +  Jl-S-+i-)  =  0.  (3-42) 

LO  L2       C2 

There  are  two  separate  classes  of  solutions  to  (3-42) 

a.  The  first  class  has  frequencies  each  of  which  exceeds  f.  In  this 
class  the  term  in  s  is  negligible,  and  to  O(s)  we  obtain  the  Poincare  modes, 
i.e., 


2    2 
K 

2 


a2  =  f2  +  C2  (k2  +  ^-)  +  O(s),  n  =  1,2,....  (  3  -  43  ) 

L 


Equation  (3-43)    shows    that  the   high-frequency   Poincare   waves   are   also 
essentially  unaffected  by  the  small  bottom  slope. 
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b.  The  class  having  frequencies  a  =  O(s),  for  which  the  first  term 
in  (3-42)  is  negligible,  while  the  second  is  of  0(1). 

The    other    solution    leads    us    to    the    dispersion     relation     for     the 
topographic  Rossby  wave,  i.e., 

a  =  -s(i) £i -,  n  =  l,2,...,  (3-44) 

,  2     n  n        r 
k  + + — 

2  2 

L  C 

obtained  by  neglecting  the  first  term  of  (3-41). 

The  maximum  Rossby-wave  frequency  should  be  given  by  setting 


22  f2       1/0 

k  =  k  ^fUL  +  JL)172,  (3-45) 

n       v       2  2 

L  C 


for  which  we  get 


G  =  a      =-^- .  (3-46) 

max  2  A  t  2 

r  2    2       IL     1/2 

[n  it    + ] 

_2 

c 


From  equation  (3-46)  it  is  clearly  seen  that  the  Rossby-wave  frequency 
is  always  less  than  f.  Another  important  feature  of  the  Rossby  wave  is  that 
its  phase  speed  in  the  x-direction,  which  is  given  by 
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s_f 

C  =f  = ,  (3-47) 

,  2      n  71         F 

L2       C20 


is  negative  for  s  >  0,  and  positive  for  s  <  0.  Therefore,  the  wave 
propagates  (in  the  Northern  Hemisphere)  parallel  to  the  topography,  with 
the  shallowest  water  on  its  right.  Also,  for  high  wave  number,  i.e.,  small 
scale,  the  frequency  a  decreases  with  increasing  wave  number  as  shown  in 
Fig.  3.3. 

The  dynamical  fields  for  the  Rossby-wave  to  lowest  order  are  given  by 


r|  =  ri0  sin(— -H  cos(kx  -  at  +  <|>)  +  O(s),  (  3  -  48  ) 


u  =  -  f^j^o  cos(^)  cos(kx  -  at  +  <J>)  +  O(s),  (  3  -  49  ) 


v  =  -  -|kT|0  sin(-^-)  sin(kx  -  at  +  (j>)  +  O(s),  (  3  -  50  ) 


where  small  terms  of  O(s)  have  been  neglected  consistently. 
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Fig.  3.3  A  schematic  representation  of  the  dispersion  diagram 
for  Poincare,  Kelvin,  and  topographic  Rossby  waves 
in  a  channel. 
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From  equations  (3-48),   (3-49),   and  (3-50),   it  follows  that  to  lowest 
order  in  s,   and  therefore  a/f,    the  fields   of  motion   in   the  Rossby   wave 

satisfy 


u=-— ^— ,  (3  -  51  ) 

f  dy 


v=   f£, 


which  are  the  time-independent  forms  of  (3-31)  and  (3-32).  Because  of  the 
close  relationship  between  rj,  and  the  pressure  field,  equations  (3-51)  and 
(3-52)  represent  the  geostrophic  relation  for  the  horizontal  motions. 

At  this  stage,  it  is  clearly  seen  that  to  lowest  order  in  a/f,  which 
plays  the  role  of  the  Rossby  number  here,  the  velocity  fields,  though 
changing  with  time,  remain  continuously  in  geostrophic  balance  with  the 
pressure  field.  However,  the  flow  is  not  exactly  geostrophic,  for  then  the 
flow  would  be  restricted  to  travel  parallel  to  the  isobaths,  i.e.,  v  would 
vanish.  Even  though  the  velocity  fields  are  geostrophic  to  lowest  order,  it 
is  the  very  small  departures  from  geostrophy  that  give  rise  to  the  wave.  It 
is  the  small  cross-isobath  flow,  which  is  a  nongeostrophic  effect,  which 
produces  the  oscillation.  The  Rossby-wave,  whose  existence  requires  both 
s  and  f  to  be  non  zero,  is  a  low-frequency  wave  oscillation;  its  period  is 
greater  than  a  rotation  period. 
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IX^KXDRAULI£_llIMPS.-in-R£TATm£ 
.AM.B.M.QM.RQIAIIM.G.S.YSTEMS. 

A.     GENERAL 

A  strong  and  relatively  warm  wind  known  as  chinook  occurs  from  time 
to  time  in  areas  along  the  eastern  slope  of  the  Rocky  Mountains.  Similar 
phenomena  are  often  observed  in  the  Owens  valley  on  the  eastern  side  of 
the  Sierra  Nevada.  Also  cold  fronts,  approaching  the  area  of  Alps,  moving 
from  north  or  west,  many  times  undergo  severe  deformation.  This 
deformation  may  result  in  a  number  of  important  weather  events  including 
blocking  and  splitting  of  air  flow  on  the  upstream  side  or  even  triggering  of 
lee  cyclogenesis  in  the  downstream  flow. 

Tepper  (1952)  has  proposed  that  squall  lines  are  modified  hydraulic 
jumps,  but  the  idea  that  there  is  a  certain  link  between  downslope  winds 
and  hydraulic  jumps  was  proposed  by  Long  (1953).  Long  suggested  that 
the  mechanism  that  produces  the  hydraulic  jump  may  be  similar  to  one  that 
produces  strong  waves  and  downslope  winds  in  the  atmosphere.  Houghton 
and  Kasahara  (1968)  have  also  proposed  other  mesoscale  hydraulic 
analogies.  Williams  and  Hori  (1970)  observed  a  delay  in  the  formation  of 
hydraulic  jumps  in  case  the  Rossby  number  was  less  than  0. 1  in  a  rotating 
system. 

However,  it  has  been  difficult  to  confirm  this  hypothesis  about  jump 
formation,  because  there  are  significant  differences  between  the  atmosphere 
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and  the  simple  fluid  systems  used  in  the  hydraulic  theory.  One  of  the  major 
reasons  for  this  uncertainty  appears  to  be  the  vertical  propagation  of  the 
wave  energy,  which  may  occur  in  the  real  atmosphere  but  not  in  fluids 
bounded  by  a  rigid  or  a  free  surface. 

B.     HYDRAULIC  JUMPS   IN  A  NONROTATING  SYSTEM 

Let  us  first  consider  the  one-dimensional  Shallow-Water  Model.  Let  us 
also  consider  a  mean  flow  over  an  isolated  rigid  orographic  ridge  as  shown 
in  Fig.  4. 1.  The  horizontal  momentum  equation  is  given  by 

du         du  d   ,,      ,     „ 

ar+uar+g*r(h+hM>=°-  f4-1* 

The  continuity  equation  also,  is  given  by 

9h         dh     ,  9u      n 

3-+u^+h^=0,  (4-2) 

dt         dx         dx 

where  h  denotes  the  depth  of  the  fluid,  and 

hM  is  the  height  of  the  rigid  ridge,  which  is  a  function  of  x. 
At  time  t  =  0,  the  fluid  is  set  in  motion  from  rest  so  that  for  infinite  x, 
it  has    a  constant   zonal   flow    uo.    After   sufficient    time    has    elapsed,    the 
solution   in   the   neighborhood   of  the  rigid  ridge  would   be   given   by    the 
steady  state  solutions  of  equations  (4-1)  and  (4-2). 
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Fig.    4.1   A  cross  section  view  of  the  one-dimensional  Shallow 
Water  Model  with  a  mean  flow  u. 
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Following  Houghton  and  Kasahara  (1968),   the  steady  state  solutions   for 
the  variables  u  and  h  are  given  by 


u2  +  2gh  +  2ghM  =  C1,  (4-3) 


hu  =  C2,  (  4  -  4  ) 


where  Ci  and  C2  are  constants.  If  a  flow  without  a  hydraulic  jump  is 
considered,  both  Ci  and  C2  are  determined  by  the  velocity  uo,  and  the 
height  ho  of  the  flow  in  the  region  far  from  the  ridge,  so  that 


Cj  =  u0  +  2ghQ,  (  4  -  5  ) 


C2  =  h0u0.  (4-6) 

At  this  stage  we  define  the  following  dimensionless  parameters  Fo,  F, 
R,  and  U*  as 


uo 
Fn  =  ^=>  (4-7) 


0 


Js\' 


F2 

2  0 

F=y-+1,  (4-8) 
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hM 

ho 


U,=  — .  (4-10) 

uo 


If  we  eliminate  h  from  equation  (4-3)  by  using  equation  (4-4),  we  are  able 
to  obtain 


(F2  -  1)  U^  +  (R  -  F2)  U,  +  1  =  0.  (4-11) 


Assuming  that  the  mean  flow  uo  is  always  positive,   then  F  >   1.   Dividing 
equation  (4-11)  by  (F2  -  1)  we  can  get 


2 
Uj  +  (^fL)  u„  +  — —  =0.  (  4  -  12  ) 


F2-' 


1  F   -  1 


Equation  (4-12)  is  a  very  important  relation  among  the  variables  U*, 
R,  and  F.  We  can  easily  plot  the  solution  U*  to  equation  (4-12)  for  given 
values  of  F2.  For  instance  if  F2  equals  to  1.02  (Fo  =  0.2),  equation  (4-12) 
becomes 


U^  +  (50R-51)U,  +  50  =  0.  (4-13) 
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We  can  consider  equation  (4-13)  as  determining  the  values  of  R  that  are 
possible  with  a  given  value  of  U*,  as  shown  in  Fig.  4.2.  It  is  interesting 
here  to  note  the  singular  behavior  near  R  =  0.  In  case  of  F  =  1.125  (Fo  = 
0.5),  equation  (4-12)  becomes 


U*  +  (8R  -  9)  U,  +  8  =  0,  (  4  -  14 ) 


which  leads  to  Fig.  4.3.  Also,  if  F  =  3.0  (Fo  =  2.0),  equation  (4-12)  can 
be  written  in  the  following  simple  form 


U*  +  (0.5R  -  1.5)  U,  +  0.5  =  0,  (  4  -  15  ) 


which  gives  Fig.  4.4. 

It  is   very   important  here   to   note,    that  if  R   has   values    lower    than 
R-critical.  where  Rcriticai  is  given  by 


Rcnticai  =  p2  -   1-8899  (F2-   1)1/3,  (4-16) 


then  three  real  roots  exist,  as  we  clearly  see  in  Figs.  4.2,  4.3,  and  4.4.  In 
the  case  where  R  is  greater  than  RCriticai>  no  physically  meaningful  solution 
exists  (no  real  roots).  In  other  words,  in  this  particular  case  we  expect  the 
occurrence  of  a  hydraulic  jump. 


59 


u 
<n 

z 

LJ 
£L 

X 

- 

Z 

z 
cr 

uJ 

> 
0 

O 


a 

UJ 

u 

5 
c 
a 

a 
UJ 

■z 


Fig.    4.2  Parameter  R  as  a  function  of  solution  U*.  as  given  by 
equation  (4-12)  for  F2  =  1.02. 
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Fig.   4.3  As  in  Fig.  4.2,  except  for  F2  =  1. 125 
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If  we  plot  equation  (4-16),  as  shown  in  Fig.  4.5,  we  are  able  to  find  three 
distinct  domains.  In  both  domains  I  and  III,  the  parameter  R  (function  of 
x),  has  values  lower  than  Rcriticai-  1°  this  case,  a  real  solution  of  equation 
(4-16)  there  exists,  which  is  physically  meaningful.  On  the  other  hand  in 
domain  II,  the  parameter  R  is  greater  than  Rcriticai>  so  no  physical  solution 
exists. 

C.     HYDRAULIC  JUMPS   IN  A  ROTATING  SYSTEM 

Although  none  of  the  numerical  experiments  concerning  hydraulic 
jumps  (results  section),  involves  rotation,  it  is  useful  to  explore  the 
particular  effects  of  rotation  on  the  formation  of  hydraulic  jumps.  The  basic 
equations  for  a  homogeneous,  one-layer,  inviscid  fluid  (Williams  and  Hori, 
1970),  are  given  by 

3u        du        dh     „ 

8r+uar+g35r-fv=0-  (4-17) 

dv         dv 

^r-+u^+fu  =  0,  (4-18) 

dt         dx 


3h        dh     L  du     n 

+  u       +h  — =0,  (4-19) 

dt  dx         dx 


where  h  is  the  depth  of  the  fluid. 
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Fig.  4.5  Classification  of  asymptotic  flow  conditions  as  a 
function  of  the  maximum  height  of  the  topographic 
ridge,  RCritical.  and  initial  flow  speed  parameter  F. 
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We  perform  a  scale  analysis  by  expressing  the  independent  variables 
as 

t  =  T  t\  (  4  -  20  ) 

x  =  Lx'.  (4-21) 

The  dependent  variables  are  also  broken  up  and  scaled  as 

u  =  U  u\  (  4  -  22  ) 

v  =  V  v\  (  4  -  23  ) 


h  =  hm+x/yh''  <4-24> 


where  hm  represents  the  mean  depth  of  the  fluid.  Also,  the  Rossby  number 
and  the  Froude  number  are  defined  as 


Ro=n:  (4-25) 

and 


f°=7tV  (4"26) 
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In  order  to  include  the  results  obtained  by  the  previous  analysis  which 
are  valid  for  a  nonrotating  system,  we  examine  the  case  where  Fo  ^  1  and 
Fq  ^  Ro-  The  appropriate  time  scale  for  doing  this  is 


T  =  —^—,  (4-27) 


and  the  appropriate  v  scale 

V  =  ^U.  (4-28) 

R0 


The  nondimensional  equations  for  this  case  are  given  by 


du'      _     ,  du'      dh'     F0     , 


+ 


V^  +  ^-^V  =  0,  (4-29) 


9t'        °     dx'     dx'    D2 

Ro 


dv'  dv' 

—  +V  —  +a'  =  0,  (4-30) 


dh'      _  ,ah'      ,,  du         du'      n 

— r+F0(u'^-+h'^-)  +  ^— =0.  (4-31) 

at  u       dx  dx  dx 


It  is  obvious  that  hydraulic  jumps  can  be  formed  through  the  action  of  the 
nonlinear  terms  in  equations  (4-29)  and  (4-31).  Even  when  Fo  is  small 
(<<1),   they  can  produce  a  jump  in  a  nonrotating  system.    If  the   Coriolis 
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term  in  equation  (4-29)  dominates  then  a  hydraulic  jump  may  be  prevented. 
This  implies  that 


F2 

—  »  F.  (  4  -  32  ) 


The  form  of  the  curve  dividing   the  jump   region   from   the   nonjump 
region  should  be  given  by 


F0  =  AR2,  (4-33) 


where  A  is  a  constant.  The  numerical  solutions  show  that  the  range  for  A  is 
from  6.0  to  7.5,  which  appears  to  be  in  agreement  with  Houghton's 
analytical  curve  corresponding  to  A  =  6.5.  For  Fo  ~  1,  the  above  scale 
analysis  does  not  apply,  since  Fo  is  then  greater  than  Ro. 

When  both  Fo  ~  1,  and  Ro  ~  1,  all  terms  in  the  equations  (4-29),  (4- 
30),  and  (4-31)  are  of  the  same  order.  In  this  case,  jumps  are  expected  to 
form.  On  the  other  hand,  if  Ro  <<  1,  the  proper  time  scale  is  1/f,  while  the 
proper  v  scale  should  be  V  =  U.  The  nondimensional  equations  can  be  then 
rewritten  as 


du'     n  r  ,  du'       1    3h'         ,     n 
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dv'  9v' 

—  +  R0u>  —  +  u<  =  0,  (4-35) 


du'     ^...Bu'^Su1      1   9u'  ,  ,  „    „,» 

3-  +  R  [u1  j—-  +  h'  ,—  +  —  ^-r  ]  =  0.  (  4  -  36  ) 

dt  u       dx  dx       F„  dx 


If  we  now  neglect  all  terms  in  Ro,  the  resulting  equations  describe  an 
inertial  oscillation  in  u  and  v.  We  do  not  expect  a  hydraulic  jump  to  form  in 
this  case  except  perhaps  after  a  very  long  time. 
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V.   MODEL  DESCRIPTION 

A.     GENERAL 

There  are  two  possible  choices  of  increasing  resolution,  where  desired 
or  required,  concerning  a  triangular  subdivision.  We  can  use  near- 
equilateral  triangles  (Cullen,  1974b)  or  equilateral  triangles  (Hinsman, 
1975).  Both  have  the  advantage  of  almost  perfect  wave  propagation 
characteristics.  For  the  same  problem  we  can  also  use  rectangular 
subdivisions,  although  it  is  not  obvious  which  subdivision  is  most 
suitable.  However,  a  major  advantage  of  the  rectangular  subdivision 
(shown  in  Fig.  5.1)  is  that  it  allows  algorithms  to  be  developed  which  take 
full  advantage  of  vector  processors.  The  interesting  point  here  is  that  the 
Galerkin  method  has  to  utilize  efficient  numerical  techniques  to  be 
considered  as  a  viable  option  for  numerical  weather  prediction. 

There  are  several  solution  procedures  available  for  the  Galerkin 
method.  No  particular  attempt  is  made  to  optimize  computational  efficiency 
in  this  research  model.  A  direct  solver  is  employed  using  a  Gaussian 
elimination  procedure.  The  matrices  from  the  Galerkin  procedure  are 
decomposed  into  upper  and  lower  block  tri-diagonal  matrices.  A 
preprocessing,  representing  the  forward  substitution  stage,  can  be  done 
once.  Any  time  a  solution  is  desired,  a  back  substitution  has  to  be 
performed.  That  is  why  the  required  coefficients  for  the  backward  step 
must    be    stored    in    a    very    efficient    manner.    This    particular    algorithm 
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Fig.    5.1    Rectangular    uniform    subdivision    for    a    channel    in 
Cartesian  coordinates. 
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represents  a  'skyline'  solver,  referring  to  the  compact  method  of  storing 
only  those  coefficients  required.  This  method  has  both  the  desired  level  of 
accuracy  and  a  high  degree  of  computational  efficiency. 

B.     EQUATION  FORMULATION 

In  order  to  integrate  the  equations  governing  the  free-surface   height 
and  velocity  of  an  inviscid  hydrostatic  incompressible  fluid  we  can  write 


|n+u|i+v|i+(M^+|i)=o  (5-D 

dt         dx        dy  dx     dy 


du        du        9u  96 

-=r-  +  U  ^—  +VTr--fv+^-=0  (5-2) 

9t         9x         dy  dx 


9v        9v        9v  36 

+  u       +v       +fu  +  «l=so  (5-3) 

dt         dx         dy  dy 


where     6  is  the  geopotential  height, 

u    is  the  east/west  component  of  the  wind, 

v  is  the  north/south  component  of  the  wind,  and 

f  is  the  coriolis  parameter. 
We  now  assume  that  the  geopotential  height  6,  in  absence  of  motion  is 
O.  Then,  in  general 
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<t>  ( x,  y,  t )  =  d>  ( x,  y )  +  (t>*  ( x,  y,  t ),  (  5  -  4  ) 

where  O  is  the  mean,  and 

(J)'   is  a  perturbation  from  the  mean.   The  governing  equations  now 
can  be  written 

(J)'t  +  OD  +  (u(l),)x  +  (v(l),)y  =  0,  (5-5) 

ut  +  <j>'x  +  Kx  -  vQ  =  0,  (  5  -  6  ) 

vt  +  <t>'y+Ky+uQ  =  0,  (5-7) 


where 


9u     dv 

D  =  T-  +  T-5  (5-8) 

ox     ay 


K  =  i(u2  +  v2),  (5-9) 


_      dv     3u 

Q  =  ^--^-  +  f.  (5-10) 


dx     dy 


Because  of  the  rapidly  moving  gravity  waves,  the  stability  condition 
for  a  numerical  integration  normally  requires  a  much  smaller  time  step  than 
for  the  simple  advection  equation  since  01/2  >>  U.   Similar  results  may  be 
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expected  with  the  the  more  complete  equations  actually  used  in  numerical 
weather  prediction.  Since  the  gravity  waves  are  usually  relatively 
unimportant  in  large-scale  weather  forecasting,  the  small  time  step  required 
for  computational  stability  increases  the  computing  time  considerably  with 
little  or  no  compensation  by  way  of  increased  accuracy,  perhaps  even  a 
loss.  On  the  other  hand,  implicit  differencing  schemes,  which  may  have  no 
restriction  on  the  size  of  the  time  step,  have  the  serious  disadvantage  of 
requiring  the  inversion  of  a  large  matrix. 

A  semi-implicit  scheme  has  the  great  advantage  of  permitting  a 
relatively  large  time  step  without  unduly  increasing  computation  time.  In 
other  words  the  semi-implicit  scheme  slows  artificially  the  propagation 
speed  of  the  fastest  gravity  waves,  which  allows  a  much  larger  time  step 
than  the  normal  Courant-Fredrich-Lewy  (CFL)  stability  criterion  and  also 
offsets  some  of  the  extra  computational  expense  required  to  solve  the 
system  of  equations  assembled  at  each  step.  For  this  reason  it  is  now 
necessary  for  the  implementation  of  the  semi-implicit  time  discretization  to 
rewrite  our  equations  in  terms  of  a  velocity  potential  %  and  a 
streamfunction  \|/  defined  by 


u  =  Xx-Vy,  (5-11) 


v  =  Xy+¥x,  (5-12) 
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with  the  following  diagnostic  relations 


y+y=u+v  (5-13) 

^xx      ^yy         x         y,  v  ' 


\i/+\i/=V-U  (5-14) 

Txx       ~yy         x        y.  v  ' 


In  other  words,  we  can  use  vorticity  and  divergence  or  velocity 
potential  and  streamfunction  as  variables  instead  of  velocity  components. 
This  means  that  second-order  derivatives  appear  in  the  equations  and 
Poisson  equations  have  to  be  solved.  Using  linear  elements,  the  scheme 
obtained  for  the  Poisson  equation  is  very  similar  to  the  finite  difference 
scheme  and  can  be  inverted  by  the  same  technique.  Cullen  and  Hall  (1979) 
showed  that  the  accuracy  of  the  Galerkin  finite  element  method  solution 
was  more  accurate  for  the  vorticity-divergence  formulation  of  the  shallow- 
water  equations  than  for  an  increase  in  resolution  with  the  primitive 
formulation.  This  unstaggered  vorticity-divergence  together  with  staggered 
variable  formulation  gives  the  best  treatment  of  geostrophic  adjustment  for 
small-scale  features  (Williams  and  Schoenstadt,  1980). 

Following  the  vorticity-divergence  approach  and  dropping  the  primes 
for  the  rest  of  the  chapter,  the  equations  become 

<j>t  +  <>D  +  (u<D)x  +  (v<|>)   =0,  (5-15) 
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Ct  +  V2<|)  +  (uQ)x  +  (vQ)   =0,  (5-16) 


Dt  +  v2(|)  +  V2  K  -  (vQ)x  +  (uQ)   =0.  (  5  -  17  ) 


where  V  2  is  the  Laplacian  operator  and  C,  is  the  relative  vorticity  given  by 


dx     dy 


Now    we    express    the    velocity    as    the    sum    of    the    rotational    and 
irrotational  components 


V  =  V    +V    ,  (5-19) 

v       x 


where 


V    =kxV\j/,  (5-20) 


Vx  =  VX,  (5-21) 


and  then  we  are  able  to  rewrite  the  equations  using 


_       9u      dv      .-.2 

D=^+37=v*-  (5"22) 
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c4-|=vv, 


obtaining  the  following 


<j,t  +  OV2x  =  -(u<|>)x-(v<|>)  (5-24) 


(V2¥)t  =  -(uQ)x-(vQ)  (5-25) 


(V  2x)t  +  V  20  =  (vQ)x  -  (uQ)y  -  (Kx)x  -  (Ky)y  (  5  -  26  ) 


Manipulating  now  the  last  equation  (5-26)  and  taking  care  of  the 
bottom  topography  (assumed  to  be  not  a  function  of  time),  we  can  easily 
obtain 


0t  +  <*>V  \  =  -  [u  (0  -  4>B)]x  -  [v  (<{»  -  0B)1  (  5  -  27  ) 


(V2y)t  =  -(uQ)x-(vQ)  (5-28) 


V2(Xt  +  <D)  =  (vQ  -  Kx)x  -  (uQ  +  K  ),  (  5  -  29  ) 


where  <j)b  is  the  bottom  topography  defined  by  the  rigid  surface  z  =  hs  (x,y) 
not  a  function  of  time  (t). 
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We  now  define  the  domain  of  integration  to  be  a  channel  encompassed 
by  solid  north-south  walls  with  east-west  cyclic  boundary  conditions.  The 
boundary  condition  at  the  walls  should  be 

V  .  n  =  0,  (  5  -  30  ) 

where  n  is  the  outward  pointing  normal  vector. 

Here  it  is  interesting  to  note  that,  in  rewriting  the  equations  in  this 
form,  we  have  increased  their  order  in  x  and  in  y  from  first  to  second  order 
and  should  expect  that  it  may  be  necessary  to  impose  further  boundary 
conditions.  However,  since  we  have  sufficient  boundary  conditions, 
already,  any  further  specification  must  not  be  arbitrary  but  should  be  a 
consequence  of  the  previous  formulations.  Along  the  walls,  the  v 
component  of  the  velocity  has  to  be  equal  to  zero,  resulting  in 


(j)y  =  -fu.  (5-31) 


The  zonal  and  meridional  components  of  the  wind  can  be  written    now  as 
u  =  -¥   +XX,  (5-32) 


v  = 


Vx  +  Xy.  (5-33) 
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Since  the  v  component  of  the  velocity  has  to  be  zero  along  the  north 
and  south  walls  then  the  obvious  boundary  condition  should  be 


¥X  +  Xy  =  0,  (5-34) 


and  we  can  satisfy  the  above  condition  by  simply  setting 

\|/  =  0,  (5-35) 

when  solving  the  vorticity  equation,  and 


Xy  =  0,  (5-36) 


when  solving  the  divergence  equation.  As  a  matter  of  fact  this  is  an 
overspecification  but  equation  (5-34)  would  be  difficult  to  apply.  Our 
initial  conditions  of  course,  must  be  specifically  selected  to  satisfy  both 
equations  (5-35)  and  (5-36). 

As  we  mentioned  before  we  use  a  semi-implicit  time  discretization 
scheme  for  reasons  of  computational  efficiency.  Basically,  this  semi- 
implicit  scheme  is  simply  a  modified  leapfrog  scheme  giving  a  net  saving  in 
the  computational  time  required  to  make  a  forecast  for  a  given  time.  The 
way  in  which  this  is  accomplished  is  to  evaluate  certain  terms  implicitly  as 
a  mean  over  times  (t  -  At)  and  (t  +  At)  rather  than  at  time  t. 
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Following  this  approach  we  evaluate  all  the  terms  on  the  left  hand  side 
of  the  equations  as  an  average  at  times  (t  -  At)  and  (t  +  At),  and  all  the  right 
hand  side  at  time  t.  The  prognostic  equations  then  become 


4>t  +  <*>[V  \  (t  -  At)  +  V  \  (t  +  At)]  =  -  [u  «>  -  <})B)]x  -  [v  ((j>  -  (()B)]  (  5  -  37  ) 


V2(\|/t)  =  -(uQ)x-(vQ)  (5-38) 


V  2[%t  +  4>  (t  +  At)  +  <»  (t  -  At)]  =  [(vQ)  -  Kx]x  -  [(uQ)  +  Ky]y .  (  5  -  39  ) 


If  we  now  solve  equation  (5-39)  for  X  (t  +  At)  and  substitute  into  equation 
(5-37)  we  can  finally  get  the  following  set  of  equations 

VV-      ♦*      =[(vQ)-K]   -[(uQ)  +  K]   ^(t'At) 

t  2       LV    ^-s  xJx      LV    ^^  yJy  9 

O  (At)  O  (At) 


V  2Y(t  -  At)         1 

+  — — -  + ([u^-^  +  tv^-Ol  (5-40) 

At  O  At  y 


V  2  ((J)*  +  %t)  =  [(vQ)  -  Kx]x  -  [(uQ)  +  Ky]y  t  (  5  -  41  ) 


V2(\|/t)  =  -(uQ)x-(vQ)  (5-42) 


where  $*  is  given  by 
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<j>  =  [  <>  (t  -  At)  -  <>(t  +  At)  ]  /  2.  (5-43) 

At  this  stage  our  initial  system  of  three  equations  in  three  unknowns 
has  been  reduced  to  two  Poisson  equations  and  one  Helmholtz  equation  to 
be  solved  at  each  time  step.  Our  first  step  concerning  the  solution 
procedure  should  involve  solving  the  <j>  equation  (5-40)  for  a  new  value  of 
(J)*.  The  second  step  should  be  then  to  solve  equation  (5-41)  for  ((f)*  +  Xt) 
and  after  substitution  for  Xt-  At  last  we  solve  equation  (5-42)  for  \\rt.  Our 
history  variables  are  <J>,  u,  and  v  and  they  are  updated  after  each  time  step 
as 

<j)  (t  +  At)  =  2(j>*  -  <)  (t  -  At),  (  5  -  44  ) 

u  (t  +  At)  =  2  At  [(xt)x  -  (Yt)y]  +  u  (t  -  At),  (  5  -  45  ) 

v  (t  +  At)  =  2  At  [(xt)y  +  (\|/t)J  +  v  (t  -  At).  (  5  -  46  ) 

In  other  words,  numerical  integration  of  the  three  forecast  equations 
involves 

a.  solve  first  a  Helmholtz  equation  for  <j), 

b.  solve  a  Poisson  problem  for  \\r,  and  finally 

c.  solve  a  Poisson  problem  for  X- 
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The  space  discretization  consists  of  expanding  the  dependent  variables 
in  terms  of  basis  functions  defined  on  a  variable  mesh,  and  then 
orthogonalizing  the  error  to  the  basis  using  the  Galerkin  procedure 
described  in  Chapter  II.  An  appropriate  approximating  function  for  the 
rectangular  subdivision  should  be  a  bilinear  function  (f).  In  this  case  the 
forecast  set  of  equations  in  Galerkin  form  become 


fry  V  f .  -  A^L_  ]  f.  =  f  {  A  [(vQ).  f.  -  4"  (K-  f ■)]  }  f- 

J  J    J     o  (At)2  J      5x  J  J    3x      J  J         ' 

-  J  f  £  [(uQ)j  fi + £  (is  9] ) f. + J ^  [£<»«>  -  v,  fj + £v«> .  Vj  fj]  1 


j* L_(t.At)ffi+fiv2x(t-At)ff.,  (5-47) 

J  •!>  (At)  J  at 


Jv  2<5j  fJ  f>  ■ "  J[£  (uQ)i  9  f> "  K (vQ)j  fi]  fi ,  < 5 " 48 ' 

-J{£R,ri»jfj+<f>jf;k  <5-49> 
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The  integral  sign  here  means  an  area  integral  over  the  domain,  the  j 
subscript  refers  to  Einstein  summation  for  the  dependent  variables  and  the  i 
subscript  is  the  ith  nodal  equation.  Following  the  integration-by-parts 
procedure  the  final  form  of  the  forecast  equations  is  given  by 


O  (At) 


J  ( £  i<»Q>i f, + k,  <£>ji} f. + !i  {w*-  w& 


[v(«  -  Vlj  # >J  }  fi  -  J  — f^  Ct  -  At)  f.  f  +  J 1  [u/t  -  At)  (^ 


9y  J 


O(At) 


+  vj(t-At)(^)j]fi-Jf0ujfjfij  (5-50) 


I  Kf »j  to  to  +  (f >j  to  to']  =  ■  J  ^ to  1  f* 


1  [^  #¥ f.  ■ 


(5-51) 
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J{J.[(vQ).,.K.(§.]}fi.J{^[(uQ)jfj  +  Kj(|)j]}f,        (5-52) 


The  line  integral  along  the  north  and  south  walls  has  been  dropped 
from  the  vorticity  equation  (5-51),  since  the  value  of  9\|//3t  is  zero  on  the 
boundaries  (\y  =  constant).  Also,  the  line  integral  along  the  north  and  south 
walls  has  been  dropped  from  the  divergence  equation  (5-52),  because  the 
value  of  the  normal  derivative  along  the  north/south  boundaries  is  zero.  As 
we  mentioned  before  the  initial  conditions  should  also  satisfy  the  condition 
that  the  normal  derivative  of  %  along  the  north/south  boundaries  is  zero. 

C.     STABILITY  ANALYSIS 

Here  we  analyze  the  primitive  form  of  the  forecast  equations  and  a 
semi-implicit  time  scheme,  since  the  results  will  be  identical  for  the 
vorticity-divergence  form,  (Hinsman,  1983).  The  one-dimensional 
equations  with  a  mean  flow,  U,  are  given  by 


9u     9<J)       IT  du 


y-U^-fo.  (5-54) 


9<b      .  9u       TT  d<b 

-^-  +  <D__  =  -U^Z-.  (5-55) 

9t         dx  9x 
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Evaluating  the  time  derivatives  with  a  centered  time  differencing,  and 
averaging  the  other  terms  on  the  left-hand  side  between  time  levels  (t  +  At), 
and  (t  -  At),  equations  (5-53),  (5-54),  and  (5-55)  become 


u(x,t  +  At)  -  u(x,t  -  At)      1    <j)(x  +  Ax,t  +  At)  -  §(x  -  Ax,t  +  At) 


2  At  2  2  Ax 


<j)(x  +  Ax,t  -  At)  -  <)>(x  -  Ax,t  -  At)  u(x  +  Ax,t)  -  u(x  -  Ax,t) 

+  2^  '.        J  " "      L~  2Ax  ] 


+  f  v(x,t),  (  5  -  56  ) 


v(x,t  +  At)  -  v(x,t  -  At)  v(x  +  Ax,t)  -  v(x  -  Ax,t)  ,.     „. 

z-r; =-U[ -— ]-fu(x,t),  (5-57) 

2  At  2  Ax 


(J)(x,t  +  At)  -  <|)(x,t  -  At)      <j)    u(x  +  Ax,t  +  At)  -  u(x  -  Ax,t  +  At) 
2At  +Y[  2Ax 


u(x  +  Ax,t  -  At)  -  u(x  -  Ax,t  -  At)  <j)(x  +  Ax,t)  -  <j>(x  -  Ax,t) 

H ^— J  -  -  U  (- — ).         {  D  -  JO  ) 

2  Ax  2  Ax 


Assuming  now  a  function,  F,  given  as 

F(x,t)  =  F'  exp[i(kx  +  cot)],  (  5  -  59  ) 

and  substituting  into  (5-56),  (5-57),  and  (5-58),  we  can  get 
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2 

4-{-Lr-[f2+(k,)2c2O]}=0,  (5-60) 

*    At2 


where 


s  =  sin(coAt)  +  (k*)  u  At,  (6-61) 


c  =  cos(coAt),  (  5  -  62  ) 


IcAx 
(k1)  =  sin(— -).  (5-63) 

Ax 


The  roots  of  equation  (5-60)  are  given  by 
s  =  0,  (  5  -  64  ) 


s2  =  (fAt)2  +  c2(kAt)2  O.  (  5  -  65  ) 


Requiring  u)  to  be  real,   the  roots  of  (5-65)   yield   the  following   stability 
criterion 


At  <  — - .  (  5  -  66  ) 

I  —  l  +  f 
Ax 
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A.     TOPOGRAPHIC  ROSSBY  WAVE 

Let  us  first  consider  the  horizontal  momentum  equations  (3-8),  and  (3- 
9).  If  we  cross  differentiate  (3-8)  with  respect  to  y,  and  (3-9)  with  respect 
to  x,  we  obtain 

3  u   3u  3u    3  u   3v  3u    9  u   9v     9  h 

3y  3t  3y  3x    3y  3x  3y  3y    -%  2   3y    3y3x' 


3  v      3u  3v        9  v     3v  9v         9  v       „  3u  3  h  , ,   „ 

^—=r-+^—^—+U ~  +  ^-^-+  V   -,      -,      H-f^-a-g^      -.       .  (6-2) 

dx  dt      dx  dx         -n  2      dx  dy         dx  dy        dx  dx  dy 


If  h  is  eliminated  from  (6-1),  and  (6-2),  then 


3     dv     3u  d     dv     3u  3     3v     3u 

3t   3x    3y  3x    3x    3y  3y    3x    3y 


„  ,3u     3u  „     3v    3u  ,  ,3u     3v  v  ,  ,    „  v 

-f(—  +  — )-(—-  — )(—  +  —  ).  (6-3) 

dx      dy        dx     dy      dx      dy 


Using  now  the  following  definition 


?so>*=ara7'  (6"4) 


86 


where  C,  is  the  vertical  component  of  the  vorticity,  equation  (6-3)  yields 


-^-  =  3f-  +  U3^  +  v3^  =  -(C  +  f)(T-  +  3-).  (6-5) 

dt       dt  dx         dy  dx      dy 


If  we  now  use  equation  (3-15),   equation  (6-4)  can  be  written  in  the 
following  form 


«.£±<*i  (6-6) 

dt        H     dt  v  ) 


or 


±«^)=o. 


where  f  is  assumed  to  be  constant. 

At  this  stage,  in  order  to  find  the  proper  initial  conditions  for  the 
particular  case  of  the  topographic  Rossby  wave,  we  can  work  out  a 
simplified  theory  for  Rossby  waves  (Phillips,  1965). 

Let  us  use  cartesian  coordinates  and  simplify  the  H  variation,  since  H 
is  a  function  of  y,  as 

H  =  D(l-sy),  (6-8) 
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where  D  is  a  constant,  and 

s  is  the  slope  of  H  in  the  y-direction. 
Expanding  the  (£  +  f)  /  H  term  as 


C  +  f_      5  +  f         (C  +  f)(l+sy)_C  +  f+Csy+£sy  ^ 


H        D(l-sy)  D  D 


and  considering  the  term  (C,  sy)  to  be  very  small,  equation  (6-7)  becomes 


^-(C  +  fsy)  =  0.  (6-10) 

dt 


Equation  (6-10)  can  be  also  written  as 


dC  dy 

dT  +  fsdH  (6-U) 


or 


4^+fsv  =  0.  (6-12) 

dt 


Assuming    small   amplitude   motion,    and    no    mean    flow,    equation    (6-12) 
becomes 
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^L+fsv  =  0.  (6-13) 


In  the  case  of  a  small  Rossby  number,  we  can  also  use  the  following 
relationships: 


5  =  V  V  (  6  -  14  ) 


and 


v-£.  (6-15) 


where  \|/  =  P  /  (pf).  Equation  (6-13)  then  becomes 


l(V2V)  +  fs(|^)  =  0  (6-16) 


or 


-(VV)  +  P(|^)  =  0,  (6-17) 


where  (3  =  f  s. 

Assuming  a  solution  of  the  form 
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V  =  \|/(y)  exp[i^i(x  -  ct)],  (  6  -  18  ) 


substitution  into  equation  (6-18)  yields  the  following  problem  for  \|/(y) 


,2 
(—7  -  M2)  (-  i|ic)  exp[i(i(x  -  ct)]  +  (3  y  (i(i)  exp[i|i(x  -  ct)]  =0  (  6  -  19  ) 

dy 


or 


,2 

— 2  -  (n2  +  &)y=  0.  (6-20) 

dy2  c 


Nondimensionalizing  the  domain  of  integration  (as  shown  in  Fig.  6.  1), 
our  boundary  problem  for  \|/(y)  becomes 


__L.(^  +  -J.)¥i  =0,  (6-21) 


dy 


with 


V|/1(0)  =  0,  (6-22) 


V1(a)=\|/2(a),  (6-23) 
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w 


v  =  a 


0.0  <a<  1.0 


v  =  0.0 


Fig.    6.1    Schematic     representation      of     the      non-dimensional 
domain  of  integration. 
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dy  (a)     d\j/2(a) 

-7—  =  — T—  '  (6-24) 

dy  dy 


for  the  lower  part  of  the  channel, and 


— -2-  -  (nj  +  -i)  v2  =  0,  (  6  -  25  ) 

dy  c 


with 

V2(D  =  0,  (6-26) 

V2(a)  =  v1(a),  (6-27) 

d^  (a)     dy2  (a) 


dy  dy 


(  6  -  28  ) 


for  the  upper  part  of  the  channel.  Here  \i\,   and  p.  2  are  the  corresponding 
wavenumbers.  Also 

Px  -Sjf,  (6-29) 

(32  =  s2  f,  (  6  -  30  ) 
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where  S\  is  the  slope  of  H  in  the  y-direction  for  the  lower  part,  and 
S2  is  the  slope  of  H  in  the  y-direction  for  the  upper  part. 
There  are  two  distinct  cases.    In  the  first  case,   the  lower  part  of  the 
channel  controls  the  phase  speed  c.  The  solutions  to  this  case  are 

V1  =  A  sinC^y),  (6-31) 

for  the  lower  part,  and 

y^Bsinht^d-y)],  (6-32) 

for  the  upper  part.  For  the  particular  case  where  a  =  0.5 

P2  =  -  P2.  (  6  -  33  ) 


7  9         Pi 

^  =  ^  +  -L,  (6-34) 


-  (Xf  =  -  2\i\  -  X*  .  (  6  -  35  ) 

Furthermore,  the  coefficients  A,  and  B  are  to  be  determined.  In  the  second 
distinct  case,  the  upper  part  controls  the  phase  speed  c.  The  solutions  to 
this  case  are 
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y2  =  Csin[A.2(l-y)],  (6-36) 

for  the  upper  part,  and 

V|/2  =  D  sinh(>.*y),  (  6  -  37  ) 

for  the  lower  part.  For  the  particular  case  where  a  =  0.5 

P2=-Pi'  (6-38) 

-  X\  =  \i\  +  -1 ,  (  6  -  39  ) 

-(\*)2  =  -2^-X,2.  (6-40) 

Once  more,  the  coefficients  C  and  D  are  to  be  determined. 

Let  us  now  consider  the  second  case,  where  the  upper  part  controls  the 
phase  speed  c.  Boundary  conditions  (6-23),  and  (6-27)  yield 

sin[X2(l  -  a)]  C  -  [expft*  a)  -  exp(-  x\  a)]  D  =  0.  (  6  -  41  ) 

In  the  same  way,  the  boundary  conditions  (6-24)  and  (6-28)  also  yield 
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-  X2  cos[X2(l  -  a)]  C  -  Xl  [exp&j  a)  +  exp(-  \  a)]  D.  (  6  -  42  ) 

Equations  (6-41)  and  (6-42)  lead  us  to  the  following  eigenvalue  problem 

*  *  * 

A. j  [exp(^1  a)  +  exp(-  Xx  a)]  sin[A.2(l  -  a)] 

+  X2  [exp(X>1  a)  -  exp(-  Xx  a)]  cos[X2(l  -  a)]  =  0.  (  6  -  43  ) 

Using  relationship  (6-40),  equation  (6-43)  becomes 


*  *  *  *  2  2  1/2 

X{  [exp(X1  a)  +  exp(-  ^  a)]  sinft^)   -  2^2]      (1  -  a)} 


*  2  2  1/2  *  * 

+  [(XX)   -  2p.2]      [exp(X1  a)  -  exp(-  %l  a)] 


*  1  7    Ml 

cosU^)   -2^2]1/'(l-a)}=0,  (6-44) 


where 


^2=~W.  (6-45) 


Here  L  is  the  horizontal  length  scale,  and  W  represents  the  width  of  the 
channel.  For  our  case  we  choose  a  channel  30°  x  30"  longitude  by  latitude, 
which  gives 
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L  =  5653510.75  m, 

W  =  4896083.93  m, 

so  |i2  equals  to  5.4413981 

We  are  able  to  solve  (6-44)  numerically  (see  Table  IV)  using  Newton'  s 
method.  Using  Table  IV,  we  can  determine  coefficients  C,  and  D  for  any 
particular  mode.  For  instance  if 

A.*  =  9.3191, 


X2  =  5.2562, 


then  in  order  to  determine  the  relation  between  C,  and  D,  the  following  set 
of  equations  has  to  be  solved 

0.4912  C  -  105.5791  D  =  0,  (  6  -  46  ) 

4.5787  C  -  984.0783  D  =  0,  (  6  -  47  ) 

resulting  in 
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A  2 
it 

* 

5.2562 

9.3191 

11.1881 

13.5791 

17.3683 

18.9967 

23.6124 

24.8347 

29.8772 

30.8523 

Table  IV  Numerical  solutions  of  equation  (6-44),  obtained  using 
Newton'  s  method. 
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„       0.4912     „ 

D  =  l05379TC  (6"47) 


or 

C=  1.0, 

D  =  0.004652  . 

Here  C  is  chosen  to  be  1.  Finally,  the  solution  to  the  case  where  the  upper 
part  controls  c,  should  be  given  by 

\|/2  =  sin  [  5.2562  (1  -  y)  ],        0.5  <  y  <  1.0,  (6  -  48  ) 

for  the  upper  part  and 

\\f2  =  0.004652  [  exp  (  9.3191  y  )  -  exp  ( -  9.3191  y  )  ], 

0.0<y<0.5,  (6-49) 

for  the  lower  part.  The  graph  of  this  solution  is  shown  in  Fig.  6.2. 

Following  the  same  approach,  with  the  use  of  Table  IV,  we  are  able  to 
obtain 
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Fig.  6.2  Graph  of  the  solution  \j/2(y)  where  the  upper  part  of 
the  channel  controls  c,  obtained  from  equations  (6-48), 
and  (6-49). 
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V2  (upper)  =  sin  [  11.1881  (1  -y)  ],        0.5<y<1.0,  (6-50) 

\|/2  (lower)  =  -  0.000716  [  exp  (  13.5791  y  )  -  exp  ( -  13.5791  y  )  ], 

0.0<y<0.5,  (6-51) 

for  the  case  where 
\  =  13.5791, 


X2=  11.1881, 


shown  in  Fig.  6.3,  and 

\|/2  (upper)  -  sin  [  17.3683  (1  -  y)  ],       0.5  <  y  <  1.0,  (  6  -  52  ) 

\\f2  (lower)  =  0.0000506  [  exp  (  18.9967  y  )  -  exp  ( -  18.9967  y  )  ], 

0.0<y<0.5,  (6-53) 

for  the  case  where 


Xx=  18.9967, 
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Fig.    6.3    As  in  Fig.    6.2,   except  for  equations   (6-50),    and   (6- 
51). 
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X2=  17.3683, 


shown  in  Fig.  6.4. 

The  corresponding  solutions  for  the  case  where  the  lower  part  controls 
the  phase  speed  c,  are  given  below 

Xj/j  (lower)  =  sin  (  5.2562  y  ),        0.0  <  y  <  0.5,  (  6  -  54  ) 

Vj  (upper)  =  0.004652  {exp[  9.3191(1  -  y)]  -  exp[  -  9.3191(1  -  y)]}, 

0.5<y<1.0,  (6-55) 

for  the  case  where 


A.*  =  9.3191, 


\  =  5.2562, 

shown  in  Fig.  6.  5, 

^  (lower)  =  sin  (  11.1881  y),        0.0  <  y  <  0.5,  (6-56) 
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Fig.    6.4    As   in   Fig.    6.2,   except  for  equations   (6-52),    and    (6- 
53). 
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Fig.  6.5  Graph  of  the  solution  ^(y)  where  the  lower  part  of 
the  channel  controls  c,  obtained  from  equations  (6-54), 
and  (6-55). 


104 


Vj  (upper)  =  -0.000716{exp[13.5791(l  -  y)]  -  exp[  -13.5791(1  -  y)]}, 

0.5<y<1.0,       (6-57) 

for  the  case  where 


X2=  13.5791, 


Xl  =  11.1881, 

shown  in  Fig.  6.6,  and 

\\f{  (lower)  =  sin  (  17.3683  y  ),        0.0  <  y  <  0.5,  (  6  -  58  ) 

\|/1  (upper)  =  0.0000506{exp[18.9967(l  -  y)]  -  exp[  -18.9967(1  -  y)]}, 

0.5<y<1.0,       (6-59) 

for  the  case  where 


X2=  18.9967, 


X1  =  17.3683, 


shown  in  Fig.  6.7. 
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Fig.    6.6    As   in   Fig.    6.5,   except  for  equations   (6-56),    and   (6- 
57). 
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Fig.   6.7    As   in  Fig.    6.5,   except  for  equations   (6-58),    and   (6- 
59). 
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In  our  experiments  we  select  initial  conditions  in  the  form  described 
by  equations  (6-54)  and  (6-55).  That  choice  appears  to  be  an  excellent  one, 
providing  us  with  the  most  proper  modes  to  observe  the  topographic 
Rossby  wave.  It  is  desired  that  the  initial  conditions  allow  a  relative 
amount  of  control  for  the  input  parameters  as  well  as  satisfing  the  boundary 
conditions. 

The  next  logical  step  is  to  determine  the  analytic  expression  for  the 
streamfunction  \\f.  Following  closely  equations  (6-54)  and  (6-55),  the  exact 
expression  for  \|/,  is  given  by 


V  =  y  [sin(^-  y)]  [sin(-^  x)]  -  Um  (y  -  yj  +  f- ,  (  6  -  60  ) 


for  the  lower  part  (0.0  <  y  <  0.5),  and 


A  *  *  2k n 

\|/  =  0.004652  —  [exp(X  y)  -  exp(  -  X2  y)]  [sin(— —  x)] 


L0 


for  the  upper  part  (0.5  <  y  <  1.0),  where 


A  =  amplitude  of  perturbation, 

W  =  width  of  the  channel  (4896083.93m), 
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L  =  length  of  the  channel  (5653510.75m), 

n  =  wave  number, 

Um  =  mean  flow  speed, 

vmid  =  middle  point  of  the  channel, 

X\  =  5.2562,  and 

X2*  =  9.3191 

The  first  term  in  the  expressions  (6-60)  and  (6-61)  represents  the 
perturbation  part.  The  second  term  is  the  north/south  slope  necessary  to 
support  a  mean  flow  of  Um.  The  last  term  plays  the  role  of  the  mean  depth 
term.  The  geopotential  height,  <}),  is  related  geostrophically  to  the 
streamfunction,  y,  by  the  following  relationship 

0  =  foV,  (6-62) 

resulting  in 

f0A 
<J>  =  -^  sin(aiy)  sin(a2x)  -  fQ  Um(y  -  y^)  +  O,  (  6  -  63  ) 

for  the  lower  part  (0.0  <  y  <  0.5),  and 
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f0A 


<{)  =  0.004652  — —  [exp(a3y)  -  exp(-  a^)]  sin(a2x) 


fnU    (y-ymJ  +  o,  (6-64) 


0      m  v/      •'mid 


for  the  upper  part  (0.5  <  y  <  1.0),  where 


»,=•£,  (6-65) 


a2  =  ^l,  (6-66) 


* 


2l1  =  X1.  (  6  -  67  ) 

The    u,    and    v    components    of    velocity    can    be    derived    using    the 
following  geostrophic  expressions 


(  6  -  68  ) 


u  = 

i  a<t> 
"lay"' 

V  = 

i  a<i> 

suiting 

in 

(  6  -  69  ) 
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A  a 

u  (lower)  =  Um  -  — -  cos^y)  sin(a2x),        0.0  <  y  <  0.5,  (  6  -  70  ) 


2 


A  a 
v  (lower)  =  —?.  sin(a1y)  cos(a2x),        0.0  <  y  <  0.5,  (6-71) 


and 


u  (upper)  =  Um  -  0.043352  —  [exp(a3y)  +  exp(-  ^y)]  sin(a2x), 

0.5<y<1.0,  (6-72) 

Aa, 
v  (upper)  =  — —  [exp(aJy)  -  exp(-  a3y)]  cos(a2x),        0.5  <  y  <  1.0.    •  (  6  -  73  ) 


The    initial    vorticity     can     be     also     derived,     using     the     following 
relationship 


C-VV.  (6-74) 


resulting  in 


C  =  -  y  [(ax)    +  (aj)  1  sin^y)  sin(a2x),  (  6  -  75  ) 


for  the  lower  part  (0.0  <  y  <  0.5),  and 
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A 


C  =  0.004652  -  [(a3)2  -  (a^2]  [exp(a3y)  -  exp(-  a3y)]  sinC^x),  (  6  -  76  ) 


for  the  upper  part  (0.5  <  y  <  1.0). 

We  set  the  initial  divergence  equal  to  zero,  assuming  the  initial  fields 
to  be  almost  geostrophic.  The  initial  fields  of  geopotential,  u-component, 
v-component,  vorticity,  and  divergence  for  the  case  where  the  lower  part  of 
the  channel  controls  the  phase  speed,  c,  and  no  topography  effect  is 
involved  (Equations  6-60  through  6-76),    are  illustrated  in  Fig.  6.8. 

B.     HYDRAULIC  JUMPS 

The  stationary  theory  predicts  the  formation  of  jumps  in  pairs,  one 
upstream  and  one  downstream  of  the  rigid  ridge  (Houghton  and  Kasahara, 
1968).  The  classification  shown  in  Fig.  4.5,  is  strictly  valid  for  non- 
rotating  flows.  No  equivalent  theory  exists  for  flows  over  mountains  in  a 
rotating  system.  Houghton  (1969)  and  Williams  and  Hori  (1970)  consider 
the  transient  motion  of  a  shallow  water  layer  on  an  f-plane  without 
mountains  starting  from  an  initial  velocity  disturbance  of  magnitude  U  over 
a  length  L. 

For  our  case  we  consider  a  nonrotating  Shallow-Water  system  (f  =  0). 
It  is  desired  once  more,  that  the  initial  conditions  allow  a  relative  amount 
of  control  for  the  input  parameters  as  well  as  satisfing  the  boundary 
conditions.  The  forecast  model  history-carrying  variables  are  <{),  u,  and  v. 
The  analytic  expression  for  the  sreamfunction,  \\f,  is  given  by 
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Fig.  6.8  Initial  conditions  for  the  GFEM  model  with  a 
rectangular  subdivision,  and  wave  number  one. 
Contour  intervals  are  600  m2/s2  for  geopotential 
height,  0.2  m/s  for  u  and  v,  0.6  x  10"6  s_1  for 
vorticity.  Nodal  points  are  denoted  by  an  x. 


113 


\\f  =  -  U  y,        0.0  <  y  <  W, 


(  6  -  77  ) 


where  U  is  the  velocity  of  the  flow  in  the  region  far  from  the  ridge.    Also 
the  initial  geopotential  height,  $o,  is  set  equal  to  the  mean  depth,  gH. 

Finally,  the  initial  u-  and  v-  components  of  the  velocity  are  given  by 


u0  =  U, 


(6-78) 


vQ  =  0. 


(  6  -  79  ) 
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XIL_EXEEEIMEttIS..AttI2_RE£iILI3. 

Our  first  experiment  involves  bottom  topography  which  is  composed 
of  two  regions.  These  two  regions  can  provide  us  with  either  an  east-west 
oriented  ridge  or  valley.  We  consider  conditions  with  no  mean  flow,  the 
objective  being  to  examine  how  well  our  model  simulates  the  topographic 
Rossby  wave,  by  comparison  with  the  theoretical  phase  speed  values. 

Our  second  experiment  is  to  investigate  the  ability  of  the  same  model 
to  create  hydraulic  jumps  analogous  to  that  predicted  from  the  analytical 
approach  (Chapter  IV).  In  this  case  we  consider  a  mean  flow  forced  to 
pass  over  a  topographic  ridge  which  extends  north-south  across  the 
channel.  In  this  experiment  we  will  consider  several  distinct  cases 
corresponding  to  different  discrete  domains  obtained  by  the  theory  (Fig. 
4.5). 

A.     EXPERIMENT  I 

We  perform  experiment  I  using  the  GFEM  rectangular  model  described 
in  Chapter  VI.  The  basic  difference  between  the  rectangular  and  triangular 
models  is  in  the  approximating  polynomials.  The  rectangular  polynomials 
are  bilinear  while  the  triangular  polynomials  are  linear.  Many  integrals 
require  evaluation  during  the  integration  process.  We  could  use  numerical 
quadrature,  however,  a  more  efficient  method  is  available  through  the  use 
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of  natural  coordinates.  We  can  accomplish  quadrature  with  no  error  by 
formula  with  this  method.  A  description  of  the  natural  coordinate  method  is 
given  in  Appendix  A  for  the  rectangular  discretization. 

The  diagram  for  the  grids  for  the  rectangular  subdivision  is  shown  in 
Fig.  5.1.  The  domain  of  integration  is  5,653.5  km  in  the  x-direction  and 
4,896.1  km  in  the  y-direction.  The  model  has  12  increments  in  the  x-  and 
y-  directions,  which  gives  the  model  156  degrees  of  freedom.  The  Ax  is 
471.1  km,  and  the  Ay  is  408.0  km.  The  value  of  the  Coriolis  parameter  is 
taken  to  be  0.00010284,  corresponding  to  45.0*  N  latitude. 

The  initial  conditions  for  the  case  of  the  topographic  Rossby  wave  are 
described  in  Chapter  VI.  A  small  wave  perturbation  is  added  to  the 
geopotential  field  which  includes  the  mean  height  and  the  required 
north/south  slope  in  case  of  non  zero  mean  flow.  It  consists  of  a  wave  with 
a  wavenumber  one,  confined  primarily  to  the  south  part  of  the  channel 
domain.  In  our  case,  we  choose  a  mean  depth  of  1,000  meters,  and  no 
mean  flow.  The  motion  is  confined  in  a  channel  with  cyclic  boundary 
conditions  as  shown  in  Fig.  7.1.  We  examine  small  amplitude  wave  motion 
encountering  different  slopes  of  the  bottom  topography  in  the  y-direction  as 
illustrated  in  Fig.  7.2,  and  7.3.  We  can  visualize  the  whole  setting  as  if  we 
placed  a  long  triangular  mountain  with  its  peak  centered  in  the  middle  of 
the  width  of  the  channel.  It  is  obvious  that  the  slope  of  the  lower  part  of 
the  channel,  S\,  corresponds  to  the  south  slope  of  the  mountain,  while  the 
slope  of  the  upper  part,  S2,  corresponds  to  the  north  slope.  No  matter  what 
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Fig.    7.1    Schematic  representation  of  the  domain  of  integration. 
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Fig.    7.2    Schematic  representation  of  the  bottom  topography  for 
the  case  of  triangular  mountain. 
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Fig.    7.3    As  in  Fig.  7.2,  except  for  the  case  of  triangular  valley. 


119 


the  case  is,  the  following  relationship  holds: 

sL  =  -s2,  (7-1) 

where   si    is   assumed   to    be    greater    than    or   equal    to    zero    for   most   of 
the  cases. 

B.     RESULTS   I 

We  integrate  the  forecast  equations  over  a  time  interval  of  96  hours, 
and  we  plot  the  results  every  48  hours.  Our  first  concern  is  to  examine  the 
case  of  no  topography.  In  this  particular  case  I,  the  peak  of  the  triangular 
mountain  is  zero,  and  both  slopes,  Si  and  S2,  are  equal  to  zero  as  we  can 
see  from  Table  V.  The  initial  conditions  are  almost  purely  geostrophic  as 
clearly  illustrated  in  Fig.  7.4.  The  integration  produces  forecast  fields 
almost  identical  to  the  initial  fields,  as  shown  in  Figs.  7.5,  and  7.6.  That  is 
expected,  since  no  topographic  or  beta  effect  is  involved.  If  we  now 
increase  slightly  the  magnitude  both  of  Si  and  S2,  as  shown  in  Table  V  for 
case  II,  the  48  hour  integration  does  not  show  any  significant  change  in  the 
forecast  fields,  as  we  can  see  in  Fig.  7.7.  However,  the  96  hour  integration 
yields  a  very  small  tendency  for  a  westward  shifting  valid  for  all  the 
forecast  fields,  and  for  the  lower  part  of  the  channel. 

Our  next  step,  case  III,  is  to  increase  the  peak  of  the  mountain  to 
163.  1  m.  In  both  the  48,  and  96  hour  integrations,  a  tendency  for  the  lower 
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Peak(<>) 

Peak  (z) 

Slope  (south) 

case  I 

0.0 

0.0 

0.0 

case  11 

400.0 

40.8 

0.0000167 

case  IE 

1600.0 

163.1 

0.0000666 

case  IV 

4800.0 

489.3 

0.0001999 

case  V 

-  400.0 

-40.8 

-0.0000167 

case  VI 

-  1600.0 

-  163.1 

-  0.0000666 

case  VII 

-  4800.0 

-  489.3 

-0.0001999 

Table  V      Peak  (in  geopotential  meters  and  in  meters),   and  south 
slope  values  for  cases  I  through  VII. 
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Fig.  7.4  Initial  conditions  for  experiment  I  (cases  I  through 
VII).  Contour  intervals  are  600  m2/s2  for  geopotential 
height,  and  0.2  m/s  for  u  and  v.  Nodal  points  are 
denoted  by  an  x. 
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Fig.   7.5    As    in    Fig.    7.4,    except    for   case    I,    and    a    48    hour 
integration. 
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Fig.   7.6    As  in  Fig.  7.5,  except  for  a  96  hour  integration. 
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Fig.   7.7    As  in  Fig.  7.5,  except  for  case  II. 
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Fig.   7.8    As  in  Fig.  7.6,  except  for  case  II. 
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part  forecast  fields  for  a  westward  shifting  is  clearly  noticed,  as  well  as  a 
small  tendency  for  eastward  shifting  corresponding  to  the  upper  part  fields 
(Figs.  7.9,  and  7.10).  In  case  IV,  which  is  our  last  case  of  positive  south 
slope  Si,  we  increase  the  peak  of  the  mountain  to  489.3  m.  This  value 
corresponds  to  almost  the  half  of  the  mean  depth  value  considered  for  our 
shallow  water  approximations.  In  both  the  48,  and  96  hour  integration  it  is 
very  evident  the  significant  westward  sifting  of  the  forecast  fields  for  the 
the  lower  part  of  the  channel,  and  the  eastward  shifting  of  the  same  fields, 
characterizing  the  upper  part,  as  we  can  see  in  Figs.  7.11,  and  7.12.  The 
results  obtained  here,  will  be  compared  quantitatively  with  the 
corresponding  analytical  values. 

At  this  stage  we  wish  to  examine  cases  V,  VI,  and  VII,  all  having 
negative  south  slopes.  The  above  mentioned  cases  correspond  to  cases  II, 
III,  and  IV,  but  with  exactly  opposite  sign  slopes.  We  can  visualize  the 
whole  setting  here  as  if  we  placed  reverse  triangular  mountains  (valleys)  of 
different  heights  with  their  lowest  points  centered  in  the  middle  of  the 
width  of  the  channel,  and  their  peak  values  always  to  be  given  by  negative 
values.  The  48  hour  integration  results  for  case  V  do  not  show  any 
significant  change  in  the  forecast  fields,  as  we  can  see  in  Fig.  7.13,  while 
the  96  hour  integration  yield  a  very  small  tendency  for  a  eastward  shifting 
valid  for  all  the  forecast  fields,  and  for  the  lower  part  of  the  channel  (Fig. 
7.14).  In  case  VI,  the  valley  is  163.1  m  deep  in  the  center.  In  both  the 
48,  and  96  hour  integrations,  a  tendency  in  the  lower  part  forecast  fields 
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Fig.  7.9  As  in  Fig.  7.5,  except  for  case  III. 
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Fig.   7.  10    As  in  Fig.  7.6,  except  for  case  III. 
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Fig.   7.11    As  in  Fig.  7.5,  except  for  case  IV. 
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Fig.   7.  12    As  in  Fig.  7.6,  except  for  case  IV. 


131 


GEOPOTENTIRL  HEIGHT 


STRERMLINES 


TM 


*     x     n 


X  X  XI 


X  X  *)       X         X  X 


X-HXIS   irETCTEl 


Fig.   7.  13    As  in  Fig.  7.5,  except  for  case  V. 
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Fig.   7.  14    As  in  Fig.  7.6,  except  for  case  V. 
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for  a  eastward  shifting  is  evident,  as  well  as  a  small  tendency  for  westward 
shifting  corresponding  to  the  upper  part  fields  (Figs.  7. 15,  and  7.  16).  The 
same  tendency  becomes  more  evident  in  both  the  48,  and  96  hour 
integration,  i.  e.,  the  significant  westward  sifting  of  the  forecast  fields  for 
the  the  lower  part  of  the  channel,  and  the  eastward  shifting  of  the 
corresponding  upper  fields,  as  it  is  clearly  shown  in  Figs.  7.  17,  and  7. 18. 
At  this  point,  we  strongly  believe  that  there  is  a  certain  link  between 
the  presence  of  topography  and  the  observed  shifting  in  all  the  forecast 
fields,  because  no  shifting  at  all  is  been  observed  in  the  flat  case  of  no 
topography.  That  specific  link  has  to  be  the  topographic  Rossby  wave 
whose  existence  requires  the  topographic  effect  in  a  rotating  system  (with 
constant  f  ).  The  reason  is  that  the  Rossby  wave  in  general,  can  exist  only 
in  the  presence  of  an  ambient  potential-vorticity  gradient.  In  case  I,  of 
course,  no  potential-vorticity  gradient  is  present,  and  that  is  why  we  do  not 
observe  any  sign  of  the  Rossby  wave.  In  case  II,  III,  and  IV,  the  positive 
y-direction  (shown  in  Fig.  7.19),  is  also  the  direction  of  increasing 
ambient  potential-vorticity.  If  we  consider  a  fluid  column  initially  at  rest, 
but  later  to  be  displaced  in  the  positive  y-direction,  then  in  order  to 
conserve  its  total  potential-vorticity,  the  wave  potential-vorticity  £o  -  Frjo, 
has  to  be  decreased.  This  will  balance  the  excess  of  the  ambient  potential 
vorticity  in  its  new  place.  There  are  two  ways  for  this  to  be  accomplished. 

a.  The  fluid  will  have  a  tendency  to  be  squeezed,    since  its  new 
position  appears  to  be  shallower  than  the  previous  one.  This  will  induce  the 
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Fig.  7. 15  As  in  Fig.  7.5,  except  for  case  VI. 
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Fig.  7.  16  As  in  Fig.  7.6,  except  for  case  VI. 
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Fig.  7.  17  As  in  Fig.  7.5,  except  for  case  VII. 
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Fig.   7.  18    As  in  Fig.  7.6,  except  for  case  VII. 
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Fig.  7.  19  The  required  ambient  potentiai-vorticity  gradient.  A 
clue  for  the  physical  explanation  of  the  topographic 
Rossby  wave  oscillation. 
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production  of  negative  vorticity  due  to  the  vortex-tube  compression. 

b.  The  act  of  the  squeezing  can  not  be  complete,  since  the  upper 
surface  is  not  restricted  or  bounded.  In  this  case  the  column  will  have  a 
tendency  to  ride  at  least  partially  up,  the  slope  resulting  in  a  greater  value 
of  rjo  in  its  new  position  than  its  neighbors. 

It  is  important  here  to  note  that  both  effects,  C,o  <  0,  and  rjo  >  0,  act  to 
reduce  the  quantity  £o  -  Frjo-  Both  effects  also,  will  result  in  a  clockwise 
circulation  in  the  fluid  around  the  column.  The  clockwise  circulation  in  the 
fluid  column  C,  will  force  the  adjacent  column  to  its  right,  R,  into  deeper 
fluid,  and  the  adjacent  column  to  the  left,  L,  to  be  squeezed  into  shallower 
fluid,  as  shown  in  Fig.  7.  20.  The  column  R,  will  become  the  center  of  a 
counter-clockwise  circulation,  while  column  L,  will  become  the  center  of  a 
clockwise  circulation.  Both  contributions  from  columns  L,  and  R,  will 
force  the  return  of  column  C,  toward  its  original  position.  This  will  result 
in  an  overshoot  due  to  the  column  C  inertia,  and  the  oscillation  will 
continue.  This  is  a  very  simplified  view  of  the  phenomenon  but  it  clearly 
shows  a  very  important  aspect.  The  strength  of  the  restoring  mechanism 
depends  on  the  vigor  of  the  circulation  induced  on  neighboring  fluid 
columns  by  the  displaced  column.  Following  the  same  approach,  in  cases 
V,  VI,  and  VII,  the  positive  y-direction  is  the  direction  of  decreasing 
ambient  potential-vorticity,  and  everything  described  above  applies  in  the 
opposite  sense. 
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Fig.  7.20  The  position  of  the  three-point  vortices  L,  C,  and  R 
.  at  three  successive  times.  Initially  collinear  and 
positioned  along  an  isobath,  C  is  displaced 
upwards,  producing  velocities  at  L  and  R  which 
move  them  as  shown.  The  vorticity  induced  on  L 
and  R  produces  a  velocity  at  C  with  a  tendency  to 
restore  it  to  its  original  position. 
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From  a  purely  theoretical  point  of  view,  if  we  recall  equation  (3-47)  we  can 
easily  find  out  that  for  positive  values  of  slope  s  the  phase  speed  in  the  x- 
direction  is  always  negative  or,  referring  to  our  case,  westwards.  For  the 
opposite  case  of  negative  values  of  s,  the  propagation  of  the  Rossby  wave 
in  the  x-direction  is  always  positive  or  eastwards.  Also,  from  the  same 
equation  (3-47),  it  is  obvious  that  for  increasing  values  of  s  (in 
magnitude),  the  propagation  phase  speed  also  increases.  In  other  words,  if 
we  keep  increasing  the  peak  of  our  triangular  mountain  we  should  expect 
the  presence  of  the  topographic  Rossby  wave  to  become  more  and  more 
evident.  In  order  to  examine  the  actual  phase  speed  in  more  detail,  we  can 
Fourier  analyze  the  v-component  field,  and  obtain  the  wave  component 
phase  speed  every  24  hours.  The  observed  phase  speed  is  then  given  by 


L  (♦,  -  <D2) 

CF  =  — -,  (7-2) 

2k  (t  - 1  ) 


where  L  is  the  length  of  the  channel.  These  phase  speeds,  averaged  over  96 
hours  for  cases  II  through  VII,  are  given  in  Table  VI,  and  in  Figs.  7.21, 
and  7.22. 

The   analytic   phase   speed   given    by    (3.41)    can    be   rewritten    in    the 
following  form: 
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Peak(z) 

C  (estimated) 

C  (corr.  fs) 

C  (corr.  fs+H) 

C  (observed) 

case  n 

40.8 

-0.72 

-0.49 

-0.49 

-0.49 

case  IH 

163.1 

-2.87 

-  1.97 

-  1.98 

-1.99 

case  IV 

489.3 

-8.61 

-6.14 

-7.12 

-7.44 

case  V 

-40.8 

0.72 

0.49 

0.49 

0.49 

case  VI 

-  163.1 

2.87 

1.97 

1.85 

1.74 

case  VTI 

-  489.3 

8.61 

6.14 

5.22 

4.89 

Table  VI  Estimated,  estimated  corrected  by  the  free  surface  term, 
estimated  corrected  by  both  the  free  surface  and  H 
terms,  and  observed  phase  speed  values  (in  m/s),  for 
cases  II  through  VII. 
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Fig.  7.21  Comparison  of  the  observed  phase  speed  values, 
with  estimated,  estimated  corrected  by  the  free 
surface  term,  and  estimated  corrected  by  both  the 
free  surface  and  H  terms,  phase  speed  values  (in 
m/s)  for  cases  II  through  IV  (scatter  diagram). 
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Fig.    7.22      As  in  Fig.  7.21,  except  for  cases  V  through  VII. 
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C  =  — — ,  (7-3) 

2*2 

~^~ 
where  h\i  is  the  mountain   height,    and   the   term  [(|ii)2   +   (X-i)2]    is    non- 
dimensional.  The  analytic  phase  speeds  resulting  from  (7-3)  are  also  given 
in  Table  VI. 

In  general,  the  analytic  phase  speeds  are  all  larger  in  magnitude  rhan 
the  corresponding  model  phase  speeds,  obtained  from  (7-2).  The  reason 
why  this  happens  is  that  the  analytic  theory  developed  here  is  based  on  a 
rigid  upper  lid,  but  at  the  same  time,  the  GFEM  model  used  for  our 
forecasts  has  an  upper  free  surface.  If  we  wish  to  include  the  upper  free 
surface  effect  in  equation  (7-3),  the  term  [(fo)2  /  gH]  must  be  added  to  the 
denominator,  resulting  in 


-fo 

1 

hM 

H 

w 

2 

2 

*? 

f2 
i      ° 

2 

w 

gH 

C  =  — : — ,  '     (7-4) 


where  H  represents  the  mean  depth  value. 
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The  value  of  the  new  added  term  is  about  the  half  of  the  denominator 
value,  so  that  it  will  reduce  the  analytic  phase  speed  value  by  around  30%, 
which  would  bring  C  and  Cp,  into  more  agreement  for  most  of  the  cases. 
However,  the  highest  mountain  peak  case,  case  IV,  or  the  lowest  depth 
valley  case,  case  VII,  can  not  be  satisfied  by  only  considering  this  change. 
The  above  mentioned  agreement  could  be  improved  for  case  IV,  by  using  a 
smaller  H  value  into  (7-4),  since  the  true  average  depth  in  this  case  is  less 
than  1  km.  Using  the  same  arguments  for  case  VII,  the  improvement  of  the 
desired  agreement  could  be  accomplished  by  using  a  larger  H  value  into  (7- 
4),  since  the  true  average  depth  in  this  last  case  is  greater  than  1  km. 

Table  VI  and  Figs.  7.21,  7.22,  7.23,  and  7.24,  compare  the  observed 
phase  speeds  with  the  various  theoretical  estimates.  As  expected,  each  new 
correction  improves  the  agreement  with  the  model  phase  speed.  The  phase 
speed  for  the  shallower  topography  are  extremely  accurate.  The  reason  for 
the  lack  of  agreement  for  larger  H  values  is  due  to  the  uncertainty  for  what 
the  correct  value  of  H  is  to  use  in  equation  (7-4).  Therefore,  for  higher 
values  of  the  bottom  topography  the  error  is  larger  because  using  the 
correct  value  of  H  in  equation  (7-4)  is  most  important.  Clearly,  experiment 
I  shows  that  our  numerical  model  is  able  to  handle  the  topographic  Rossby 
wave  case  extremely  well. 
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Fig.    7.23       As  in  Fig.  7.21,  except  for  bar  representation. 
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Fig.   7.24       As  in  Fig.  7.22,  except  for  bar  representation. 
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C.     EXPERIMENT  II 

We  perform  experiment  II  using  the  same  GFEM  rectangular  model,  as 
in  experiment  I,  the  only  difference  is  in  the  resolution  of  the  model.  The 
model  now  has  24  increments  in  the  x-  and  y-  directions  which  gives  the 
model  576  degrees  of  freedom.  The  domain  of  integration  is  727.6  km  in 
the  x-direction,  and  630.  1  km  in  the  y-direction.  The  value  of  Coriolis 
parameter  is  taken  to  be  zero  for  all  of  the  cases,  corresponding  to  a 
nonrotating  system. 

The  initial  conditions  for  experiment  II  are  described  in  Chapter  VI. 
The  boundary  conditions  in  x  are  chosen  to  be  periodic  once  more,  that  is 

a  (0)  =  u  (L),  (  7  -  5  ) 

(J)  (0)  =  4>  (L).  (7-6) 

These  boundaries  at  x  =  0,  and  x  =  L,  are  placed  sufficiently  far  from  the 
ridge  so  that  the  desired  asymptotic  conditions  are  well  established  in  the 
vicinity  of  the  ridge  before  wave  motions  are  able  to  be  fed  back  into  this 
region  by  the  periodic  boundary  conditions.  The  height  profile  of  the 
orographic  ridge  is  given  by 


{ 


HMsin2^}    for  0^x^z' 


hM=    1  (7-7) 

0    elsewhere, 
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where  Z  is  the  width  and  Hm  the  height  of  the  mountain,  as  shown  in  Fig. 
7.25. 

D.     RESULTS  II 

We  integrate  the  forecast  equations  over  a  maximum  time  interval  of 
50.6  minutes  and  we  plot  the  results  at  t  =  5.6,  16.9,  28. 1,  39.4,  and  50.6 
minutes.  The  five  distinct  cases  we  run  (I  through  V),  are  summarized  in 
Tables  VII  and  VIII.  Cases  III  and  IV  are  expected  to  produce  a  hydraulic 
jump  because  they  lie  in  domain  II  (see  Fig.  4.5).  In  each  case  the 
equations  are  integrated  from  an  initial  state  where  u  and  h  are  both 
uniform.  The  u-field  for  case  I  is  shown  in  Fig.  7.26.  It  clearly  shows  the 
rapid  development  of  a  speed  maximum  over  the  center  area  of  the  mountain 
and  slightly  on  the  lee  side.  A  secondary  speed  maximum  also  forms  over 
the  ridge  area  and  it  moves  upstream  with  time.  The  <J>- field,  shown  in  Fig. 

7.27,  indicates  the  earlier  development  of  low  pressure-field  over  the  lee 
side  of  the  ridge.  The  high  pressure-field  over  the  east  part  of  this  pattern 
has  an  obvious  tendency  to  move  upstream  with  time.  Since  the  main 
feature  of  the  flow,  the  speed  maximum  over  the  center  area  of  the 
mountain,  is  quite  stable  we  regard  this  as  a  no-jump  case  in  agreement 
with  the  theory  presented  before  in  chapter  IV.  The  u-field  for  case  II,  a 
case  with  slightly  higher  mountain,  and  stronger  mean  flow,  shown  in  Fig. 

7.28,  is  generally  similar  to  case  I  except  that  the  perturbation  now 
represents  a  larger  fraction  of  the  mean  flow.  The  same  arguments  appear 
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Fig.  7.25  Schematic  representation  of  the  position  of  bottom 
topography  along  x-axis,  valid  for  each  node  per 
horizontal  row,  for  the  hydraulic  jump  case 
(experiment  II). 
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Fo 

R 

F 

H(m) 

U(m/s) 

D 

case  I 

0.2 

0.1 

1.010 

1000.0 

19.80 

I 

case  II 

0.4 

0.2 

1.040 

1000.0 

39.60 

I 

case  IE 

0.3 

0.7 

1.220 

1000.0 

29.69 

n 

caselV 

0.8 

0.2 

1.149 

500.0 

56.00 

n 

case  V 

1.4 

0.05 

1.407 

400.0 

87.65 

in 

Table  VII  Froude  number  (Fo),  maximum  height  of  the  ridge  (R), 
parameter  F,  mean  depth  (H),  mean  flow  (U),  and 
domain  (D),  for  cases  I  through  V. 
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F 

R 

D 

Jump  case  ? 

case  I 

1.010 

0.10 

I 

no 

case  II 

1.040 

0.20 

I 

no 

case  HI 

1.220 

0.70 

II 

yes 

case  IV 

1.149 

0.20 

II 

yes 

case  V 

1.407 

0.05 

III 

no 

Table    VIII 


Parameter  F,  maximum  height  of  the  ridge  (R), 
domain  (D),  and  classification  of  the  asymptotic 
flow  conditions,  for  cases  I  through  V. 
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to  be  valid  for  the  <J)-field  also,  shown  in  Fig.  7.29.  We  also  regard  case  II 
as  a  no-jump  case  in  agreement  with  the  theory. 

The  u-field  for  case  III  is  shown  in  Fig.  7.30.  In  this  case  the  wind 
maximum  on  the  lee  side  of  the  mountain,  continues  to  grow  with  the  time. 
At  the  same  time  the  0-f ield,  shown  in  Fig.  7.31,  indicates  that  the  low 
pressure-field  centered  on  the  lee  side,  changes  rapidly  to  nearly  zero.  We 
regard  case  III  as  a  jump  case  because  it  is  not  approaching  steady  state. 
The  resolution  of  the  model  is  too  poor  to  allow  a  detailed  description  of 
the  small-  scale  jump  zone.  The  theory  also  indicates  that  this  is  a  jump 
case.  Fig.  7.32  contains  the  u-field  for  case  IV.  The  <f)-field  for  the  same 
case  is  given  by  Fig.  7.33.  The  behavior  of  case  IV  is  similar  to  case  III, 
so  that,  this  is  also  a  jump  case  in  agreement  with  the  theory. 

The  fields  of  u  and  0  for  case  V,  are  shown  in  Figs.  7.34  and  7.35 
respectively.  These  field  patterns  are  similar  to  those  of  cases  I  and  II,  the 
only  difference  is  a  downstream  shifting  for  both  the  u-  and  <})-  fields.  The 
growth  in  the  amplitude  of  the  curves  appears  to  be  stabilized,  so  that  we 
regard  case  V  as  a  no-jump  case  in  agreement  once  more  with  the  theory. 

In  all  the  investigated  cases  (I  through  V),  jump  formation  is  indicated 
by  both  u  and  <J)  amplitudes  which  continue  to  increase  with  time.  In  order 
to  study  the  behavior  of  each  one  of  the  jump  cases  in  more  detail,  better 
resolution  will  be  required  as  well  as  a  larger  domain  to  reduce  the 
boundary  effects. 
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In  this  thesis  we  have  investigated  how  well  a  particular  shallow-water 
finite  element  prediction  model  handles  surface  topography.  In  both 
experiments  the  flow  is  confined  within  an  east-west  channel  with  periodic 
boundary  conditions. 

In  the  first  experiment  the  bottom  is  composed  of  two  regions  of 
constant  and  opposite  north-south  slope  with  no  east-west  variation,  so  that 
the  bottom  is  either  an  east-west  ridge  or  an  east-west  valley.  In  order  to 
have  the  proper  initial  conditions,  the  analytic  Rossby  wave  solutions  are 
derived.  These  are  obtained  by  solving  the  linearized  quasi-geostrophic 
equations  with  the  free  surface  assumed  to  be  rigid.  Each  solution  has  a 
sinusoidal  variation  with  y  over  one  bottom  slope  and  exponential  decay 
over  the  other  slope.  With  the  simple  Rossby  formula  the  direction  of 
propagation  is  determined  by  the  bottom  slope  in  the  region  which  has  the 
largest  wave  amplitude.  These  solutions  are  similar  to  the  trench  wave 
solutions  obtained  by  Mysak  et  al  (1979),  who  used  piecewise  exponential 
bottom  profiles.  The  numerical  solutions  with  the  initial  conditions  given 
by  the  linear  solutions  produced  smoothly  propagating  solutions.  The  phase 
speeds  were  very  accurate  when  the  Rossby  formula  was  corrected  for  free- 
surface  effects  and  mean  depth. 

In  the  second  experiment  a  north-south  ridge  was  placed  across  the 
channel  with  sine-squared  east-west  variation.  The  Coriolis  parameter  was 
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set  to  zero  and  the  initial  u  and  <J>  were  constant.  The  theory  of  Houghton 
and  Kasahara  (1968)  for  the  formation  of  hydraulic  jumps  was  reviewed. 
The  equations  were  integrated  for  five  initial  conditions.  In  each  case  a 
speed  maximum  occurred  over  or  just  downstream  from  the  ridge  and  low 
heights  were  found  on  the  lee  side  of  the  ridge.  For  three  of  the  cases  the 
solutions  reached  an  approximate  steady  state.  These  agreed  with  the  theory 
which  predicted  no  jumps.  The  two  other  cases  lead  to  increasing  winds 
and  decreasing  heights  with  no  steady  state.  These  were  jump  cases 
according  to  the  theory.  In  these  two  last  cases  the  model  resolution  was 
inadequate  to  simulate  the  formation  of  the  hydraulic  jumps  in  detail. 

Finally,  the  finite  element  model  performed  well  for  two  very  different 
topographic  effects.  Further  testing  is  required  on  the  jump  cases  with 
higher  resolution.  Furthermore,  the  effect  of  the  semi-implicit 
discretization  on  the  hydraulic  jumps  should  be  determined. 
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APPENDIX 
NUMERICAL  QUADRATURE 

A  very  fast  and  efficient  method  is  required  for  the  evaluation  of  the 
integrals  obtained  by  the  Galerkin  approach.  For  the  rectangular 
subdivision,  integration  formulas  are  based  on  an  orthogonal  axis 
transformation.  In  this  case,  the  integrals  to  be  evaluated  contain  either 
products  of  the  basis  functions,  products  of  derivatives  of  basis  functions, 
or  a  mixture  of  both. 

An  orthogonal  axis  transformation  will  allow  quadrature  formulas  for 
rectangles.  We  are  able  to  transform  the  rectangle  shown  in  Fig.  A.  1  by 
using  the  following 


C"*.  (A-l) 

a 


n=^.  (A-2) 


where  the  values  of  £,   and  rj    at  each   corner  are   shown   in   parentheses. 
Using  fj  as  a  basis  function,  we  can  express  f,  as 


(1+CiOd+TliTl) 

f  •  = ; ! •  (  A  -  3  ) 
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Fig.    A.  1  Orthogonal  axis  transformation  for  rectangular 

integration  formulas. 
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Derivatives  of  the  basis  functon  are  given  by 

df.      1    df. 

^  =  --^  (A-4) 

ox      a   ^r 

and 

df.      1    df. 

ay    b  ori 

The  interaction  coefficients  can  now  be  determined  for  a  derivative  or 
straight  inner  product.  The  straight  inner  product  is  given  by 

1   i 


Cij  =  |J  fifjdxdy  =  abJ|ffjdCd7i 


-l  -l 


=  ~  Jd+CiOd  +CfC)dCjd  +T]iTl)(l  +TljT|) 


dTl 


~  (2  +  |-CC)(2  +  |-Tl.Tl).  (A-6) 


Also,  the  mixed  derivative  is  given  by 
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r  r  df.  df.  k   r  r  3f-  df  ■ 

d-=     I  3J-3J-dxdy  =  -   [  |-JL-id?dn 

1J   J  J  ^  ^  *  J  J  ac  ac 


i  i 


b 
16a 


JCiCjdCj(l+Tl.Tl)(l+TljTl)dn 


■  1  -1 


16a         *  J  3     *    J 


(A-7) 
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